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LOW  REYNOLDS  NUMBER  CALCULATION  OF  TURBULENT  CHANNEL 
FLOW:  A  GENERAL  DISCUSSION 


1.  Introduction 

Direct  numerical  simulations  (DNS)  of  turbulent  shear  flows  have  recently  provided 
researchers  with  a  new  tool  that  has  been  used  to  study  the  dynamics  of  turbulence 
in  great  detail.  In  some  cases  these  calculations  have  provided  information  that  would 
be  difficult  or  impossible  to  obtain  from  experiments  and  in  other  cases  they  have 
provided  new  results  that  were  latter  confirmed  by  experiments.  At  the  Naval  Research 
Laboratory,  we  have  recognized  that  DNS  may  provide  new  insights  into  problems  that 
are  of  particular  interest  to  the  Navy.  Problems  to  which  DNS  may  be  applied  in 
the  future  include  investigations  of  the  sources  of  turbulent  flow  noise,  the  physics  of 
polymer-induced  drag  reduction,  and  turbulence-free  surface  interactions.  However, 
in  order  to  embark  on  a  detailed  study  of  any  of  these  complex  problems,  it  is  first 
necessary  to  establish  confidence  in  our  ability  to  perform  meaningful  simulations  of 
turbulence  in  relatively  simple  cases. 

The  channel  geometry  has  been  chosen  by  many  investigators  as  a  starting  point  for 
DNS  wall-bounded  shear  flow  calculations  because  of  both  its  geometrical  simplicity  and 
the  ease  with  which  it  can  be  handled  numerically.  However,  in  spite  of  its  simplicity, 
intense  investigations  of  the  high  resolution  calculations  of  Kim.  Moin  and  Moser(  1987) 
(KMM)  .  for  example,  have  yet  to  unequivocably  unravel  the  spatial  and  temporal 
evolution  of  the  coherent  vortical  structures  which  experimental  evidence  suggests  have 
their  origin  near  the  wall.  Therefore,  the  complex  physics  of  even  this  simple  wall- 
bounded  shear  flow,  it  seems,  warrants  continued  intensive  investigation  in  the  hope 
that  such  studies  may  reveal  more  details  of  fundoinont.nl  wall  turbulence  phenomena 
of  which,  perhaps  the  most  important  is  the  nature  of  the  Reynolds  stress  producing 
events.  We  have  therefore  undertaken,  in  the  last  few  years,  an  effort  in  DNS  which 
para  llels.  in  some  respects,  the  efforts  of  other  groups  for  the  purpose  of  eventually 
applying  these  methods  to  attack  the  more  complex  problems  mentioned  above.  In 
this  report  we  review  the  progress  that  has  recently  boon  made  at  NRL  in  moderate 
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resolution,  low  Reynolds  number  DA'S  channel  how  calculations.  In  the  future  we  will 
report  on  progress  that  has  been  made  in  applying  DXS  to  more  complicated  geometries 
and  on  the  results  of  higher  resolution  calculations. 

In  Section  2  we  review  the  numerical  methods  that  have  been  used  to  perform  DNS 
calculations  of  turbulence.  In  particular,  we  discuss  the  difficulties  associated  with  oper¬ 
ator  splitting  methods  and  the  methods  that  have  been  devised,  such  as  Green  function 
methods,  for  overcoming  operator  splitting  errors.  We  also  propose  a  modification  of 
the  original  Green  function  method  devised  by  Marcus(19S4)  which  we  will  analyze 
and  test  in  future  investigations.  We  also  discuss  methods  that  have  been  devised  to 
solve  the  equations  of  motion  without  resort  to  splitting.  We  conclude  Section  2  with 
a  review  of  the  computer  codes  that  we  have  used  and  note  the  algorithms  upon  which 
each  is  based.  We  have  created  numerous  datasets  which  cover  differing  length  and 
time  scales  which  can  be  used  for  distinctly  different  purposes.  The  attributes  of  each 
dataset  is  listed  at  the  end  of  Section  2. 

In  Section  3  we  discuss  one  and  two  point  statistics  from  several  datasets  and  also 
present  some  preliminary  conditional  sampling  results.  We  find  that  in  spite  of  the  low 
resolution  of  our  calculations,  our  results  compare  well  in  most  respects  to  the  higher 
resolution  calculations  of  IvMM  and  also  with  experiments.  We  note,  however,  that 
our  computations  give  somewhat  lower  intensities  and  somewhat  longer  streamwise 
length  scales  than  the  KMM  results.  We  speculate  on  several  possible  reasons  for  these 
differences. 

\\  e  also  include,  in  Appendix  A.  a  discussion  of  t  he  effects  of  aliasing  errors  on  the 
long-time  behavior  of  the  calculations.  We  present  evidence  that  the  numerical  solu¬ 
tions  to  the  Xavier-Stokes  equations  may  be  severely  and  adversely  affected  by  aliasing 
errors.  The  evidence  presented  indicates  that  if  these  errors  are  not  removed,  the  solu¬ 
tion  reaches  a  chaotic  steady  state  which  is  remarkably  similar  to  a  weak  transitional 
turbulent  flow.  That  is.  we  have  found  that  aliasing  errors  have  a  damping  effect  on 
the  turbulence  which  has  yet  to  be  understood. 


2.  A  Review  of  Numerical  Methods 


Here  we  review  the  principal  numerical  methods  that  have  been  employed  in  the  di¬ 
rect  numerical  simulation  of  turbulence  by  numerous  investigators  over  the  past  decade. 
We  place  special  emphasis  on  a  review  of  the  methods  that  have  been  used  in  the  sim¬ 
ulation  of  wall  bounded  turbulence.  Orszag  and  Patterson  (1973)  performed  the  first 
simulation  of  homogenous  isotropic  turbulence  which  required  323  grid  points.  A  spec¬ 
tral  Galerkin  approach  was  employed  in  which  trigonometeric  functions  were  used.  The 
next  major  direct  simulation  was  the  turbulent  channel  flow  calculation  of  Orszag  and 
Kells(lOSO)  and  Orszag  and  Patera(  19S3).  In  these  simulations  one  homogenous  coor¬ 
dinate  using  trigonometric  functions  was  replaced  by  C'hebyshev  polynomials.  In  the 
first  paper,  the  problem  of  transition  to  turbulence  in  a  channel  geometry  was  studied, 
and  in  the  second,  the  study  was  extended  to  turbulent  flow.  At  the  same  time  there 
was  work  being  done  at  NASA  Ames  (Kim  and  Moin(19S2))  and  other  institutions  on 
large  eddy  simulations  (LES).  LES  involves  many  of  the  same  numerical  techniques,  but 
includes  a  term  in  the  Navicr-Stokes  equations  which  models  the  unresolved  turbulent 
length  scales. 

The  next  major  step  in  direct  simulations  was  taken  by  Marcus  (19S4)  and  Schu¬ 
mann  (19S5).  This  improvement  involves  the  use  of  Green  functions  or  so-called  in¬ 
fluence  matrix  methods  to  improve  the  implementation  of  the  boundary  conditions. 
This  method  improves  the  convergence  properties  of  the  earlier  methods  and  will  be 
discussed  in  detail  shortly.  Marcus  applied  this  method  to  the  Taylor-Couctte  problem 
and  satisfactorily  predicted  wave  speeds.  Schumann  applied  his  version  of  this  method 
to  turbulent  channel  flow. 

There  are  difficulties  in  the  above  methods  related  to  satisfying  the  incompressibil¬ 
ity  condition  which  cannot  easily  be  overcome.  An  alternative  is  to  rewrite  the  Navier- 
Stokes  equations  as  a  fourth-order  equation  for  the  normal  velocity  and  a  second-order 
equation  for  the  normal  vorticity.  In  this  form,  the  equations  are  more  difficult  to  solve, 
but  inrompressibilfy  is  satisfied  implicitly. 


2.1  Governing  Equations  and  the  Channel  Flow  Geometry 

We  are  particularly  interested  in  studying  turbulent  flow  between  two  rigid  parallel 
walls.  This  choice  represents  a  compromise  between  scientific  interest  and  numerical 
effort.  The  streamwise  coordinate  is  x  or  j-j ,  the  wall  normal  coordinate  is  y  or  x^ 
and  the  spanwise  coordinate  is  z  or  £3.  The  computational  domain  used  throughout 
this  work,  in  terms  of  channel  half-widths,  is  2  in  the  wall  normal  direction  and  5 
in  both  the  streamwise  and  spanwise  directions.  All  calculations  are  performed  with 
C’hebyshev  polynomials  in  the  wall  normal  direction  and  Fourier  series  in  the  streamwise 
and  spanwise  directions. 

The  governing  equations  for  this  problem  are  the  Navier-Stokes  equations  in  the 
so-called  rotational  form: 

<9u  u  •  u  1 

'at  =  U  X  V(i>+  T~)  +  u  +  f>  (2.1) 

and  the  continuity  equation, 

V  ■  u  =  0.  (2.2) 

In  these  equations,  p  +  ^  is  the  dynamic  pressure  head,  Re  is  the  Reynolds  number 
based  on  channel  half-width  and  friction  velocity,  and  f  is  a  body  force. 

It  is  possible  to  eliminate  the  pressure  from  the  above  equations.  This  results  in 
the  fourth  order  equation: 
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for  the  normal  velocity.  In  addition  there  is  a  second  order  equation  for  the  normal 
vorticitv: 

r)  1 

(2.4) 
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Fr'2  ~  A-  +  r?  u’’2' 


with 
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In  this  formulation  the  flow  conserves  mass  implicitly. 

For  all  the  schemes,  no-slip  boundary  conditions  are  imposed  at  the  walls.  From 
the  continuity  equation  it  follows  that  an  additional  boundary  condition  at  the  walls, 


du  2 

d.r> 


=  0  on  xo  =  ±1, 


can  be  derived.  This  boundary  condition  is  needed  when  solving  the  fourth  order 
equation  for  the  normal  velocity.  In  the  streamwise  and  spanwise  direction  the  use  of 
trigonometric  functions  requires  periodic  boundary  conditions. 

There  are  several  ways  to  non-dimensionalize  the  Navier-Stokes  equations.  The  two 
methods  used  currently  are  inner  scaling  and  outer  scaling  whose  definitions  are  given 
in  Table  2.1.  When  describing  turbulent  physics  close  to  a  wall  the  only  parameters 
available  that  can  Ire  used  to  form  length  and  time  scales  are  the  viscosity  and  the 
wall  shear  stress.  For  inner  variable  scaling  we  define  the  friction  velocity,  u*,  a  length 
scale,  /*,  and  a  friction  or  wall  Reynolds  number,  R* .  This  scaling  is  the  appropriate 
one  when  discussing  the  flow  behavior  in  the  viscous  sublayer,  the  buffer  layer  and  to 
some  extent  the  logarithmic  layer. 

\\  hen  descibing  the  physics  far  from  the  wall,  outer  scaling  is  used.  In  turbulent 
channel  flow  the  outer  scales  are  usually  taken  to  be  the  mean  centerline  velocity,  Uci, 
and  the  channel  half-width.  L>.  Outer  scaling  is  appropriate  when  discussing  the  large 
scale  properties  of  the  flow.  In  the  current  state  of  direct  simulation,  the  core  region 
of  the  flow  is  small  and  this  scaling  h  not  as  effective  as  inner  scaling.  In  turbulent 
transition  studies  outer  scaling,  however,  is  more  appropriate. 


Table  2.1 

Non-dimensionalizations 


Variable 

Inner  Scaling 

Outer  Scaling 

Velocity 

"*= y-'S  i . . 

I'd 

Length 

/*  = 

L> 

Viscosity 

V 

V 

Reynolds 

R*  rzz  SlLz 

Rc  = 

V 

Number 

When  inner  scaling  is  used,  the  driving  force  can  be  related  to  the  mean  wall  shear- 
stress  (and  the  friction  velocity)  in  the  sense  that  the  wall  shear  stiess  will  converge 
asymptotically  to  or  rather  oscillate  about  a  mean  value.  The  wall  shear  stress  and  the 
friction  velocity  can  be  related  to  the  driving  force  in  the  following  manner: 


Ou 

dx> 


<n\ 

(lx> 


where  P()  is  the  static  pressure.  A  large  scab'  velocity  is  not  specified  since  it  is  not 
needed.  When  outer  scaling  is  used  this  relationship  must  hold,  but  there  is  the  ad 
ditional  need  to  specify  a  large  scale  velocity.  Since  there  is  no  unique  relationship 
between  the  mean  turbulent  centerline  velocity  and  the  wall  shear  stress,  the  most  nat¬ 
ural  velocity  to  choose  is  the  initial  mean  centerline  velocity.  Since  our  main  interest 
is  fully  developed  turbulence  and  not  turbulent  traur  ifion.  all  further  discurNun:  will 
be  in  the  context,  of  inner  scaling. 
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2.2  Time  Stepping  Methods 


The  incompressible  Xavier -Stokes  equations  form  a  non-linear  coupled  set  of  partial 
differential  equations.  There  are  two  major  difficulties  with  these  equations  when  trying 
to  select  a  time-stepping  scheme:  the  non-linearity  of  the  convective  term  and  the 
enforcement  of  the  incompressibility  condition.  These  two  problems  are  handled  in 
very  different  ways.  The  extreme  complexity  of  the  non-linear  term  makes  treating  it 
implicitly  very  unattractive  so  that  it  is  generally  handled  with  an  explicit  or  so-called 
semi-implicit  numerical  method. 

The  second  problem,  enforcing  incompressibility,  is  equivalent  to  solving  the  time 
dependent  Stokes  problem.  The  Stokes  problem  is  the  Xavier- Stokes  problem  with  the 
incompressibiltv  equation  but  with  no  convective  term.  The  statement  of  this  problem 
is  given  as  follows: 


and 


3u 

~dt 


vP+7LvV, 


(2.5) 


V  •  u  =  0 


(2.0) 


u  =  0  on  .r2  =  ±1.  (2.7) 

There  are  two  methods  for  dealing  with  this  problem:  the  split  schemes  and  the 
unsplit  scheme.,.  I:i  the  split  schemes  the  coupled  linear  Stokes  equations  are  separated 
into  a  pressure  step  and  a  viscous  step.  The  incompressibility  condition  is  enforced 
during  the  pressure  step.  In  the  unsplit  schemes  a  fourth  order  equation  for  the  vertical 
velocity  is  formed.  It  satisfies  the  incompressibility  condition  implicitly.  Xote  that 
even  the  unsplit  scheme  requires  the  splitting  off  of  the  non-linear  terms.  The  ‘split.'  or 
’ uiispl it  character  of  a  scheme  in  the  context  of  this  paper  refers  only  to  the  algorithm 
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ust’d  ro  solve  the  Stokes  problem.  A  variety  of  split  schemes  are  used  in  the  computer 
codes  for  turbulence  simulation  now  operational  at  XRL  and  a  code  using  the  unsplit 
scheme  is  now  under  development. 

2.2.1  Time  Split  Methods 

The  motivation  for  using  the  split  scheme  is  numerical  simplicity  and  efficiency.  It 
involves  reducing  a  complicated  set  of  coupled  equations  into  a  simpler  set.  Unfortu¬ 
nately.  the  simplification  occurs  at  the  expense  of  an  improper  implementation  of  the 
boundary  conditions.  In  the  Stokes  problem  described  above,  the  pressure  field  exists 
to  satisfy  the  incompressibility  condition.  Each  of  the  time-dependent  equations  is  a 
second  order  equation  with  two  boundary  conditions:  the  zero  velocity  condition  at 
i‘2  =  ±1.  They  form  a  coupled,  bur  closed  set  of  equations.  We  do  not  wish  to  solve 
this  coupled  set  of  equations. 

The  Stokes  problem  can  be  simplified  by  splitting  it  in  two  parts  as  follows: 

The  incompressibility  stop  is. 


Ou 

Of 


=  -  V/-. 


(2.8) 


V  u  =  0. 


ami  the  viscous  step  is 
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(2.9) 


(2.10) 


u  -  0  on  j- 1  =  il. 


(2.11) 


I  lie  splitting  scheme  introduces  distinct  time  stepping  errors  ( ( )rs/ng.  Israeli  and  Dev- 
ille  198 1  )  which  dominate  flaw  of  the  viscid  or  inviscid  problems.  Therefore  there  is 


no  advantage  in  specifying  a  high  order  time  stepping  scheme.  In  addition,  the  steep 
spatial  gradients  near  the  walls  require  a  fine  mesh  and  these  terms  must  by  handled  im¬ 
plicitly  in  order  to  avoid  a  strong  time  step  constraint.  Typically,  an  Adams-Bashforth 
scheme  is  used  for  the  nonlinear  terms  and  an  implicit  Euler  method  for  the  linear 
terms.  A  Crank-Xicholson  scheme  may  also  be  used. 

We  note  that  boundary  conditions  have  not  yet  been  specified  for  the  incompress¬ 
ibility  step.  The  only  proper  boundary  condition  is  the  no-slip  condition  imposed  after 
the  pressure  step  (Marcus.  19S4).  However,  a  boundary  condition  must  be  specified  at 
the  intermediate  step.  In  this  step  there  are  four  equations  for  four  unknowns  and  they 
can  be  solved  for  a  single  variable  which  is  either  pressure  or  normal  velocity.  Solving 
for  the  pressure  field  results  in  a  Neumann  problem  while  solving  for  the  normal  ve¬ 
locity  leads  to  a  Dirichlet  problem.  There  is  no  difficulty  in  specifying  the  boundary 
conditions  for  the  linear  viscous  step. 

A.  Dirichlet  Methods 


The  original  Orszag-Kells(19SO)  algorithm  for  simulating  turbulence  between  flat 
plates  uses  a  normal  velocity  formulation  to  satisfy  the  incompressibility  condition.  The 
semi-implicit  Adams-Bashforth  Crank-Xicholson  scheme,  which  will  be  described  later, 
is  used  to  treat  the  non-linear  term.  A  backward  Euler  scheme  is  used  in  the  pressure 
and  viscous  steps.  Due  to  a  strong  memory  constraint,  an  explicit  Adams-Bashforth 
scheme  is  used  in  all  the  algorithms  discussed  in  this  work,  although  the  treatment  of 
the  non-linear  term  is  still  under  study.  When  the  non-linear  terms  are  included,  the 
three  steps  in  the  algorithm  can  be  written  as  follows: 


(a  )fhe  non-linear  step: 


u  -  11" 

At 


(2.12) 


where 


F  =  u  x  u  +  f 
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(b)  the  pressure  step: 


„,2  ,  O2  O2  0,0.  0  . 

'  M2  —  (  TTT  +  TTT  ) 11 i  ~  W - (  W - «1  -I"  7t - M3  ) 


and  (c)tlie  viscous  step: 


Or2  Ox2'"*  dx2'dxi 
u  2  —  0  on  j'2  =  ±1. 


U,,+  1  -  u  _  J_V2U»  +  , 


3x: 


(2.13) 


u 


A/ 

»+\ 


Rc 

0  on  .r  2  —  ±1. 


(2.14) 


Note  that  only  during  the  intermediate  pressure/continuity  step  is  the  continuity  equa¬ 
tion  satisfied,  while  only  at  the  end  of  the  viscous  step  are  the  boundary  conditions 
properly  satisfied.  This  inconsistency  can  be  seen  by  examining  the  normal  momentum 
equation  of  the  continuity  step: 


f/o  =  f/o  -  At— rin+1. 
dx  2 


(2.15) 


where  II  =  p  +  which  implies  that  at  the  walls 

fnn+l  =  ~  on  x  2  =  ±1. 
OX)  At 


(2.16) 


This  is  not  consistent  with  the  Xavier-Stokes  equations  for  the  normal  momentum 
evaluated  a t  the  walls  and  has  the  effect  of  introducing  spurious  pressure  forces  at  the 
walls. 

B.  Neumann  Methods 


A  boundary  condition  lor  the  pressure  field  can  be  developed  based  on  the  as¬ 
sumption  that  a t  high  Reynolds  numbers,  the  inviseid  boundary  conditions  are  good 
approximations.  From  the  Xavier-Stokes  equation  in  the  wall-normal  direction  the 
following  boundary  condition  can  be  developed: 


o 

0X2 


n  = 


1  c V 


n<  Ox 2 


on  x  >  =  ±1. 


(2.17) 
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which  in  the  inviscicl  limit  becomes 


d 

dx, 


n  =  o  on  x>  =  ±1. 


The  comiilete  algorithm  is  given  by: 
(a)  the  non-linear  step: 


(b)  the  pressure  step: 


and  (c)  the  viscous  step: 


^  _  — P"  —  1 

At  ~  2  2 


u 


V-IT+1  =  V  •  — 
At 


an 


n+ 1 


a.i'2 

U'«+i-a  i 


-  0  011  .T‘2  =  ±1. 


2,  n+1 


V-u 


At  i?c 
u',+1  =  0  on  ,r2  =  ±1. 

From  (2.20).  the  pressure  boundary  condition  implies 


it  2  =  u  2  on  x  2  =  ±1. 


(2.18) 


(2.19) 

(2.20) 


(2.21) 


(2.22) 


We  find,  therefore,  that  there  is  a  fictitious  flow  through  the  boundary  during  the 
incompressibility  step.  One  way  to  avoid  this  difficulty  is  by  using  a  Green  function 
approach. 

C.  Green  Function  Methods 

The  difficult}-  with  the  Orszag- Kells  algorithm  is  the  lack  of  proper  boundary 
conditions  at  the  intermediate  pressure  step.  In  an  attempt  to  minimize  the  errors 
induced  by  these  difficulties,  Marcus  (1984)  developed  a  Green  function  method  to  add 
a  correction  to  the  calculated  pressure  field.  He  solved  the  incompressibility  step  for  the 
pressure  and  made  the  observation  that  the  pressure  step  is  a  linear  non-homogeneous 
problem  with  non-homogeneous  boundary  conditions.  The  solution  can  then  by  divided 
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into  three  parts:  one  which  satisfies  the  inhomogeneous  equation  with  homogeneous 
boundary  conditions  and  two  other  solutions  which  satisfy  the  homogeneous  equation 
with  non- homogeneous  boundary  conditions.  Let 

nn+1  =  n,,+1 +n;,+1  (2.23) 

where  fl'i+1  is  the  solution  to 

AtV2nn  +  1  =  V  •  u,  (2.24) 


subject  to 


fr  +  1  =  0  on  .r,  =  ±1, 


and 

v2nc  =  0 


(2.25) 

(2.26) 


subject  to  inhomogenous  boundary  conditions  at  the  upper  and  lower  boundaries.  The 
boundary  conditions  for  the  pressure  correction  can  be  found  by  substituting  (2.23) 
into  (2.20)  and  evaluating  them  at  the  walls  subject  to 


.i?+1 


0  on  ,r?  =  ±1. 


(2.27) 


This  leads  to  the  following  condition: 


A/~n 

OX) 


<i  +  i 


ii i  +  t~-V2u"+1,  011  x2  =  ±1. 
Itc 


or 


=  U)  +  ^V2(/"+1  -  A/^n"+\  on  x>  =  ±1. 
cJ.7‘2  iu  Ud'2 


(2.2S) 


(2.29) 


Note  that  ii.';  1  is  not  known  during  the  incompressibility  step.  However,  Marcus 
asserts  that  this  boundary  condition  is  true  if  and  only  if: 


V  -  u',+  1  =  0. X-2  =  ±1. 


(2.30) 
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To  satisfy  these  boundary  conditions,  the  Laplace  equation  for  the  pressure  corrections 
is  solved  for  arbitrary  Dirichlet  boundary  conditions.  The  solution  can  be  written  as 
(following  Marcus'  notation): 


n"  +  1(.T2,  kXi ,  kx 3)  =  a”+1(kXl ,  kXa)\i(xi,l'xl,kXa)  +  ci%+1(kXl ,  kX3)x2(x2,kXl,kX3) 

(2.31) 

where  \'i  and  \2  are  solutions  to  the  Laplace  equation  with  the  following  boundary 


conditions: 


f  1,  *2  -  l; 
Xl  [0,  .X‘2  =  1, 


(2.32) 


cm*.* 

/  0,  .r2  =  1; 

X2  (  1,  .r2  =  -1. 

Note  that  they  are  independent  of  the  fluid  flow  and  need  to  be  calculated  only  once. 
Then  let 


where 


u  =  v  +  llc, 


c,  =  u  -  At  vn,,+1, 


(2.33) 


(2.34) 


L  =  -  At  vn"+1. 


The  final  viscous  step  is 


(2.35) 


(2.36) 


where 


and 


";,+1  =  fl  - 
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The  operator  1  —  ^V2  is  inverted  using  the  proper  Dirichlet  boundary  conditions. 
The  remaining  undetermined  coefficients  and  a 2  +  1  are  found  by  applying  the 

boundary  condition  given  by  (2.27). 

D.  A  Proposed  Alternate  Green  Function  Method 

The  Green  function  formulation  developed  by  Marcus  is  only  a  correction  to  the 
splitting  algorithm  and  still  suffers  from  one  of  the  main  difficulties  associated  with 
splitting  schemes.  Since  the  Navier-Stokes  equations  consist  of  three  second  order 
operators  (three  momentum  equations)  they  form  a  sixth-order  system  and  therefore 
require  six  boundary  conditions  for  their  solution.  In  the  time  splitting  algorithms 
there  are  usually  four  second-order  operators  with  a  total  of  eight  boundary  conditions 
which  must  be  solved.  The  Marcus  correction  improves  the  splitting  algorithm,  but 
does  not  correct  this  problem. 

The  following  algorithm  is  an  attempt  to  solve  this  problem.  The  proposed  algo¬ 
rithm  requires  six  boundary  conditions  and  only  three  Poisson  equation  inversions.  It 
has  not  been  implemented  to  the  best  of  our  knowledge  nor  has  the  algorithm  been  an¬ 
alyzed  for  stability  or  accuracy.  The  algorithm  is  similar  to  that  of  Marcus  in  its  use  of 
Green  functions.  It  is  included  here  since  it  may  be  the  basis  of  transition  calculations 
using  non- Newtonian  constitutive  equations.  In  the  next  section  an  algorithm  will  be 
presented  which  is  better  than  the  split  algorithms  but  it  should  be  recognized  that,  for 
complex  non-linear  constitutive  equations  in  non-Newtonian  flows  a  spli t.  scheme  may 
be  more  manageable. 

The  objective  is  to  solve  for  the  normal  velocity  in  the  continuity  and  viscous  steps 
using  both  zero  flow  and  zero  flux  boundary  conditions  applied  after  the  viscous  step. 
This  can  be  accomplished  through  the  use  of  Green  functions  which  couple  these  two 
stops  through  the  boundary  conditions.  The  Green  functions  must  be  evaluated  and 
the  accompanying  coefficients  must  be  determined  to  satisfy  the  boundary  conditions. 
The  updated  streaniwise  and  spanwise  velocities  can  then  be  computed  at  the  conti¬ 
nuity  step.  No  boundary  conditions  will  be  imposed  at  this  step.  The  remaining  two 
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equations  in  the  viscous  step  for  the  spanwise  and  streamwise  velocity  can  be  combined 
to  form  an  equation  for  the  normal  vorticity.  Solving  for  the  normal  vorticity  requires 
only  one  Poisson  inversion  and  two  boundary  conditions.  The  boundary  conditions  can 
be  derived  from  the  remaining  no-slip  boundary  conditions.  The  application  of  two 
sets  of  boundary  conditions  in  the  normal  velocity  problem  is  not  improper  since  two 
sets  of  second  order  problems  are  being  solved.  The  result  is  the  same  order  as  the 
original  problem  and  there  is  no  need  to  invent  boundary  conditions  for  the  pressure 
field.  By  solving  the  problem  in  this  fashion,  a  boundary  condition  for  the  pressure 
field  is  never  specified.  The  use  of  a  no  flux  boundary  condition  is  consistent  with 
enforcing  continuity  at  the  wall.  The  accuracy  and  stability  of  this  scheme  has  not  yet 
been  determined.  The  Green  function  code  currently  being  used  will  be  modified  to 
test  this  algorithm. 

2.2.2  Time  Unsplit  Methods 

The  unsplit  schemes  are  preferable  in  that  they  satisfy  continuity  implicitly  and 
the  normal  velocity  and  normal  vorticity  are  uncoupled  if  the  non-linear  terms  are 
neglected.  The  fourth-  order  problem  is  solved  using  several  auxiliary  problems  and 
the  second-order  problem  for  vorticity  is  solved  directly.  This  formulation  was  first 
discussed  by  Orszag  and  Patera  (19S3)  and  later  implemented  in  a  simpler  form  in 
the  large-scale  direct  simulations  by  Kim,  Moin  and  Moser  (1987).  It  is  the  latter 
formulation  that  will  be  used. 

Using  the  Crank-Nicholson  scheme  for  the  linear  terms  and  the  Adams-Bashforth 
scheme  for  the  non-linear  terms,  (2.3)  and  (2.4)  can  be  rewritten  as 


(1  - 


At 

Wt 


V2)V 


V,+1  - 


At 


At 


=  a  +  )V ""  +  T(3^  -  N:r l). 


(2.37) 


and 


At 


(1  -  ^V2K+1  =  (1  +  ^-V2K  +  =^(3iV"  -  AT"-1 ), 


At 


At 


2  Re 


2  Rt 


UJ  2 


U>2 


(2.38) 
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where  the  known  terms  have  been  moved  to  the  the  right  hand  side.  These  equations 
are  subject  to  the  boundary  conditions: 


uo  =  ~  ~  =  u.’)  -  0  on  x  >  =  ±1. 

'  dx2 


(2.39) 


The  vorticity  problem  is  an  inhomogeneous  Poisson  problem  with  homogeneous  bound¬ 
ary  conditions.  It  can  be  solved  relatively  easily.  A  Chcbyshev-Tau  method  will  be 
used.  The  vertical  velocity  problem  can  be  solved  if  it  is  broken  into  a  homogeneous 
problem  and  two  inhomogenous  problems: 


+  a+i/"+1  +a~u'!_+i. 


(2.40) 


which  satisfv  the  conditions 


u„  —  u+  =  u_  =  0  oxr  X‘2  =  ±1, 


(2.41) 


- — («„  -f  o  +  «  +  +  a~u_)  =  0  on  x2  —  ±1. 

OX  2 

To  solve  these  problems,  introduce  the  auxiliary  variable  (j  such  that 


(2.42) 


C+1  =  Y2t/2n+l. 


(2.43) 


Then  the  particular  problem  is 


11  -  mv‘K"’+' -,l  +  +  f,3A’" 


-  A*"'1 


(2.44) 


subject  to 


v2//;;+1  =c;+1. 


=  o  on  .i-j  =  ±i. 


(2.45) 


(2.40) 


C  =  0  on  r  ,  =  ±1. 


1G 


Then  the  particular  problem  for  satisfying  the  upper  boundary  condition  is 


(1  oRt  V2)C++1  =  °, 

(2.48) 

and 

V2  u"+1  =  C++1, 

(2.49) 

subject  to 

1  =  0  on  x2  =  ±1, 

(2.50) 

with 

,-n  +  i  f  1,  if  x2  =  1; 

^  \  0,  if  x2  =  —  1. 

(2.51) 

Then  the  particular  problem  for  satisfying  the  lower  boundary  condition  is 

(1-^v2)C-+1=0’ 

(2.52) 

and 

v2u”+1  =  c"+1, 

(2.53) 

subject  to 

u'L+1  =  0  on  x2  =  ±1, 

(2.54) 

„+1  f  0.  if  x2  =  1; 

^  "\l,  if  x2  =  — 1. 

(2.55) 

Given  the  solutions  to  the  auxiliary  problems,  (2.42)  can  be  used  to  find  the  unknown 
coefficients  a+  and  a~. 


2.3  The  Non-Linear  Term 

Treatment  of  the  non-linear  term  in  the  Navier-Stokes  equation  has  two  difficulties. 
The  first  is  the  evaluation  of  its  spatial  character.  Originally.  Orszag  and  Ivells(  1 9S0 ) 
suggested  the  use  of  the  so-called  rotational  form  of  the  Navier-Stokes  equations  in 
which  the  non-linear  term  is  written  as  ux  w.  The  objective  was  to  recast  the  non-linear 
term  in  a  quasi-conservative  form  as  suggested  by  Fornberg(  1973 ).  If  the  calculation 
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is  not  dealiased  this  form  may  not  be  quasi-conservative.  As  discussed  in  Appendix  A, 
the  use  of  this  term  with  and  without  the  aliasing  terms  leads  to  significantly  different 
results.  Horiuti  ( 1987)  and  Zang(lOSS)  have  studied  this  problem  numerically  and  both 
suggest  that  if  the  calculation  is  not  de-aliased  then  there  arc  bcttci  forms  for  the  non¬ 
linear  term  than  the  rotational  form.  All  the  simulations  currently  being  preformed  at 
XRL  are  de-alaised. 

As  discussed  earlier,  there  are  two  different  time  stepping  schemes  used  for  the 
non-linear  term.  If  the  non-linear  portion  of  the  time  split  Xavier- Stokes  equations  is 
written  as 


du 

dt 


=  U  X  u>, 


(2.50) 


then  straight  forward  application  of  the  Adams-Bashfortli  scheme  is 


^  ^  _  — Fn _ ip'1-1 

At  _  2  2 


(2.57) 


where  F  is  the  non-linear  term.  In  order  to  reduce  the  convective  instability  of  this  step 
Orszag  and  Kells  (19S0)  suggested  rewriting  the  non-linear  step  by  adding  U(x 2)57^ 
to  both  sides  of  (2. 50).  The  left  hand  side  is  treated  using  the  Crank-Nicolson  implicit 
scheme  and  the  right  hand  side  by  the  Adams-Bashforth  scheme.  The  resulting  time 
stepping  equations  are: 


u  —  u" 

A/ 


|F"  -  ~Fn_1  -  in,,)(^ 

2  2  2  Ch'i 


,<9a"  du"-1 

2 - | - 

d.i'i  d.n 


(2.58) 


Note  that  the  semi-implicit  scheme  requires  an  additional  six  variables.  If  the  interme¬ 
diate  results  u"  and  u"-'  are  not  used  the  order  of  accuracy  is  reduced  to  O(At).  Due 
to  the  large  memory  requirement,  this  scheme  is  not  usually  used. 


2.4  Spectral  Methods 

Here  we  give  only  a  brief  summary  of  spectral  methods  necessary  to  describe  the 
different  computer  codes  currently  at  XRL.  The  reader  is  referred  to  several  complete 
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reviews  of  spectral  methods  (Gottlieb  and  Orszag,  1978,  Gottlieb,  et  al.,  1984  and 
Canuto.  et  al.  1988)  for  more  detailed  descriptions. 

The  Navier-Stokes  equations  are  discretized  in  the  planes  parallel  to  the  walls 
using  a  Fourier  pseudo-spectral  scheme.  In  the  direction  normal  to  the  channel  walls 
Chebyshev  polynomials  are  used.  Neumann  or  Dirichlet  boundary  conditions  can  be 
imposed  at  the  wall.  In  the  flow  direction  and  in  the  spanwise  direction  periodic 
boundary  conditions  are  used  and  this  facilitates  the  use  of  a  trigonometric  Fourier 
series  in  these  directions.  The  velocity  field  is  represented  as 

A/-1  A'-i  P 

»<*•'>=  ZEE  u(m,n,p.t)exp(i(am  +  /3n))Tp{x  2),  (2.59) 

m=  —  A/  n=  —  .V  />= 1 

where  o  =  -P1-,  9  =  and  i  =  \f—l- 

The  algorithms  described  above  require  at  least  four  Poisson  solutions.  For  exam¬ 
ple  in  the  Dirichlet/normal  velocity  formulation  a  Poisson  equation  is  inverted  for  the 
normal  velocity  in  the  pressure  step  and  three  Poisson  equations  are  inverted  during 
the  viscous  step  to  update  the  three  velocities.  In  the  operational  computer  codes  at 
NRL  two  different  methods  are  used  in  inverting  the  Poisson  equation.  They  are  the 
collocation  and  the  tau  methods.  In  the  collocation  method  an  interpolating  polyno¬ 
mial  is  used  to  convert  continuous  equations  in  space  to  discrete  equations.  It  can  be 
shown  that  this  is  formally  equivalent  to  taking  a  discrete  cosine  transform  of  the  data 
at  the  collocation  points.  If  C  is  the  matrix  representation  for  the  discrete  Fourier 
transform,  C-1  its  inverse  and  D  the  discrete  derivative  operator  in  matrix  form  (Got- 
tleib,  et  al.  19S4).  then  the  matrix  operator.  M,  which  takes  a  real  space  function  into 
the  derivative  of  the  real  space  function  is: 

M  =  C-‘DC 

.  and  the  second  derivative  operator  is  simply 

MM. 
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The  boundary  conditions  can  be  built  into  this  matrix  and  the  solution  to  the  Poisson 
equation  is  the  solution  to: 

M-u  =  s. 

where  u  is  the  solution  and  s  is  some  source  term.  Since  MM  is  a  full  matrix,  normal 
numerical  methods  are  not  effective.  Haidvogel  and  Zang  (197S)  describe  a  method  of 
decomposing  the  full  matrix  into  eigenvalues  and  eigenvectors.  Then  the  inversion  be¬ 
comes  a  simple  matter  of  matrix  multiplication.  After  the  decomposition,  this  method 
is  fast  and  retains  the  infinite  order  accuracy  of  spectral  methods. 

In  the  Tan  method,  the  solution  to  a  problem  with  its  boundary  conditions  is 
represented  with  n  +  k  expansion  functions.  The  coefficients  are  found  using  n  equations 
derived  from  the  governing  equations  and  k  equations  derived  from  the  k  boundary 
conditions.  The  solution  is  represented  in  Fourier- C'hebyshev  coefficients.  In  the  case 
of  the  ,v>  direction,  the  Poisson  equation  can  be  rewritten  as  a  quasi -tridiagonal  system 
that  can  be  efficiently  solved  with  a  modified  Gaussian  elimination  algorithm.  The 
boundary  conditions  are  explicitly  built  into  the  matrix  operator. 

2.5  Computer  Codes  and  Data  Sets 

There  are  four  operational  computer  codes  for  the  direct  simulation  of  turbulence 
in  a  channel  currently  in  use  at  NRL.  The  first  two  listed  in  Table  2.2  were  written 
by  Orszag  and  his  students  using  the  Orszag-Kells  algorithm  (Dirichlet  boundary  con¬ 
ditions  in  the  pressure  step).  In  the  first  program,  a  very  clever  method  is  used  to 
manage  the  memory  and  input-output,  (i/o).  Only  a  quarter  of  the  computational  do¬ 
main  is  in  computer  memory  at  any  time,  and  to  prevent  wasted  i/o  time  tfiis  time 
had  to  lie  minimized  and  carefully  overlaid  by  concurrent  epu  time.  In  the  second  and 
remaining  computer  codes,  the  data  is  always  resident  in  the  computer  memory.  The 
last  two  operational  computer  codes  were  developed  by  John  McLaughlin  (Clarkson 
Lniversity).  and  extensively  updated  by  and  obtained  from  Steven  Lyons  and  Thomas 
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Hanratty  (University  of  Illinois).  The  fourth  computer  code  uses  the  Green  function 
algorithm  developed  by  Marcus. 

Table  2.2 

Computer  Codes  For  Turbulent  Channel  Flow 


Number 

Split /I'nsplit 

Boundary  Conditions 

Memory 

Status 

Author 

1 

Split 

Diriohlet 

Out  of  Core 

Operational 

t 

2 

Split 

Dirichlet 

In  Core 

Operat  ional 

+ 

+ 

3 

Split 

Neumann 

In  Core 

Operational 

0 

4 

Split 

Neumann  with 

In  Core 

Operat  ional 

0 

Green's  Functions 

5 

Split 

Neumann  with 

In  Core 

Under 

o 

Green’s  Functions 

Development 

6 

Fourth  Order 

Dirichlet  and  Neumann 

In  Core 

I’nder 

• 

Development 

T  Orszag,  Kells  and  Patera 
j  Orszag,  Bullister  and  Pelz 
0  McLaughlin 
O  Leighton  and  Handler 
•  Handler,  Leighton.  Swean  and  Wang 

There  are  two  channel  codes  under  development  at  NRL.  The  first  is  an  adaptation 
of  code  4  using  the  alternate  Green  Function  algorithm  described  in  Section  2.2.1.d. 
The  second  and  more  important  one  is  based  on  the  fourth  order  algorithm.  Much  of 
the  work  for  both  has  been  completed,  and  both  should  be  in  the  testing  phase  by  the 
end  of  calender  year  1988. 


Currently,  there  are  four  distinct  datasets  which  exist  or  will  exist  in  the  near 
future.  Each  of  these  datasets  may  have  more  than  one  entry  which  will  depend  on 
the  resolution  in  time.  There  are  an  additional  two  datasets.  A  and  DAG.  The  former 
was  a  high  resolution  simulation  which  still  contained  aliasing  errors.  See  Handler, 
Leighton  and  Carroll(  19S7)  for  a  discussion  of  this  dataset.  Program  4  was  run  for 
a  short  time  (  lOOf*).  in  order  to  develop  data  for  determining  mean  statistics.  The 
short  run  time  will  have  an  effect  on  the  mean  statistics,  hut  is  long  enough  to  get  a 
fair  idea  of  the  flow  properties.  The  remaining  datasets  are  to  be  written  to  tape  using 
the  PDSDL'MP  routine  on  the  Cray  and  can  be  brought  up  when  needed. 

In  the  datasets  C’han'2  and  C'han3.  the  difference  in  time  scales  between  entries  was 
chosen  to  be  about  a  factor  of  ten.  The  C'hani.l.  i=  1-3  dataset  are  used  for  documenting 
the  long  time  averages,  evaluation  the  proper  orthogonal  decomposition,  and  similar 
processes  when*  statistically  independent  realizations  are  needed.  The  Chani.2  datasets 
are  for  use  in  conditional  sampling,  flow  visualization  and  for  following  the  general 
evolution  of  the  flow.  Dataset  C'hani.l  and  Chani.3  .  i  =  2.3  are  for  following  the  details 
of  temporal  evolution. 


Table  2.3 

Data  Sets  For  Turbulent  Channel  Flow 


Name 

Domain  f 

(y.z.x) 

Resolution 

(  Ny.Nz.Xx  1 

R* 

Afo 

Program 

St  at  us 

Notes 

A 

2.5.0 

05. 128.0  t 

1 05 

1 

deleted 

Highly  Aliased 

Chan  1  1 1 

2.5.5 

33,6-1.16 

125 

100 

3 

exists 

28  realizations 

C It  a 1 1 1 . 2 

2.5,5 

33.61.16 

125 

l 

3 

exists 

2800  realizations  O 

DAO 

2.5.5 

33,61.16 

125 

! 

4 

exists 

Diagnostics  only 

(  han2. 1 

2.5,10 

33.61.32 

125 

50 

1 

exists 

50  realizations 

Chan2.2 

2,5,10 

33,61.32 

125 

2 

1 

exists 

500  realizations 

Chan2.3 

2.5,10 

33.6-1.32 

125 

0.2 

1 

exists 

128  realizations 

Chan3. 1 

2.5.10 

65.61,61 

150 

50 

1 

planned 

50  realizations 

Chan3.2 

2.5,10 

65.61,61 

150 

2 

4 

planned 

500  realizations 

( 'han3.3 

2.5,10 

65,61,61 

150 

0.2 

■1 

planned 

128  rcalizat  ions 

Chant .  1 

2.5,5 

65.61,32 

125 

0.1 

4 

exists 

200  realizations 

T  Based  on  channel  half  width,  .r  >  X  .i’i  X  .i'3 . 

O  Planar  Data  Only. 

O  l  ime  Between  realizations,  all  datasets  three  dimensional  except  Chan!  . 2. 
I  A  Do  called  I  ).\  ill  text. 

;  See  Handler,  heighten  and  Carroll  (l!)s7) 
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3.  Results 

3.1  One-point  statistics 

In  this  section  we  will  compare  the  one-point  statistics  computed  from  the  the 
dealiased  (DA)  and  the  dealiased-Grcen  function  (DAG)  datasets,  whose  characteristics 
have  been  disscussed  in  Section  2.5,  with  experimental  data  and  with  the  KMM  results. 
For  any  given  component  of  velocity  (or  vorticity).  which  we  will  designate  as  v  for 
convenience,  we  compute  the  mean  value,  77.  tin*  mean-square  value,  u2.  the  skewness, 
a3,  and  the  flatness,  a4,  as  follows: 

A' 

77(.r2)  =  l/.V^T  <  a-'  >,  (3.1) 

j= t 
:Y 

ft^U'2 )  =  l/.v£  <  (1  lj-  <  uJ  >)2  >,  (3.2) 

7=1 

.V 

^(.v2)  =  1  /.Y]T  <  (a3-  <  uJ  >f  >,  (3.3) 

7=1 

and 

.v 

ID  km )  =  l/A'^  <  (a3  -  <  uj  >)4  >,  (3.4) 

7  =  1 

where  <>  signifies  averaging  over  a  horizontal  plane  and  the  index  j  identifies  each  of 
the  X  realizations  of  the  flow.  Unless  otherwise  noted,  we  use  the  friction  velocity,  a*, 
and  the  viscous  length.  /*.  as  defined  in  Section  2.1  to  make  all  variables  in  this  and 
subsequent,  sections  nondirncnsional.  We  use  the  convention  that  any  variable  with  a 
raised  asterisk  has  been  made  nondimensional  using  viscous  units. 

A  typical  turbulent  velocity  profile  should  exhibit  four  distinct  regions  which  are 
named,  in  order  of  increasing  distance  from  the  wall:  the  sublayer,  buffer  layer  .  loga¬ 
rithmic  region  .  and  the  core.  It  can  be  shown  from  theoretical  considerations,  (see  for 
example  Tennekes  and  Lmnley.  1975)  that  the  sublayer  velocity  profile  must  bo  of  the 
form  J/j*  =  .7’ 2 *  and  that  a  logarithmic  region  must  exist  for  which  it\*  =  Aln( j-2*  )  +  B. 
Hussain  and  Reynolds  (1975)  assert  that  in  fully  developed  channel  flow  .4  ==  2.44  and 
D  =  5.0  although  there  is  some  scatter  in  the  experimental  estimates  given  for  these 
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constants.  The  mean  velocity  profile  obtained  from  the  DA  dataset  is  presented  in 
Figure  1  along  with  the  profiles  from  Hussain  and  Reynolds  and  KMM.  We  find  that 
the  DA  dataset  gives  good  agreement  with  the  experimental  result  and  the  higher  res¬ 
olution  calculations.  The  logarithmic  layer  in  the  current  calculation  is  not  as  thick  (in 
viscous  units)  as  in  t lie  KMM  result.  This  is  a  direct  consequence  of  the  wall  Reynolds 
number  difference  (125  for  DA  and  DAG  and  ISO  for  KMM)  since  the  wall  to  the  cen¬ 
terline  distanc  in  viscous  units  is  equal  to  the  wall  Reynolds  number  itself.  It  does  not 
appear,  however,  that  the  logarithmic  layer  from  the  DA  dataset  is  as  clearly  distinct 
from  the  buffer  layer  as  in  the  experiments  or  the  KMM  calculations. 

The  intensity  profiles  for  the  DA  dataset  are  shown  in  Figure  2.  The  peak  in 
the  intensity  of  the  streamwise  component  of  velocity  occurs  at  a  distance  of  14.76 
for  the  DA  data  which  is  in  excellent  agreement  with  KMM  and  the  experiments  of 
Kreplin  and  Eckelmann(1979)  but  the  peak  intensity  level  is  somewhat  lower  than 
either  KMM  or  experiment.  The  DA  and  the  KMM  spanwise  and  wall-normal  velocity 
intensities  are  both  considerably  lower  than  the  experimental  results.  (KMM  reference 
the  work  of  Perry.  Lim  and  Henbest(  19S5)  who  maintain  that  cross-contamination 
in  hot  wire  measurements  may  be  the  cause  of  the  higher  wall-normal  and  spanwise 
velocity  intensities  measured  by  many  investigators.)  Farther  from  the  wall  the  DA 
results  show  slightly  lower  intensities  than  the  KMM  results.  We  also  note  that  there 
is  an  unusually  high  value  of  the  wall  normal  velocity  intensity  at  a  distance  of  2.4. 
This  is  an  artifact  of  the  application  of  a  Neumann  boundary  condition  on  the  pressure 
in  the  Orszag-Kells  algorithm  which  has  been  discussed  in  Section  2.2. l.b. 

The  mean  velocity  profile  and  turbulence  computed  from  the  DAG  dataset,  which 
was  generated  using  the  Green  function  method  discussed  in  Section  2.2.I.C.  are  pre¬ 
sented  in  Figures  3  and  4.  Tin  DAG  mean  velocity  profile  exhibits  a  more  clearly 
defined  logarithmic  layer  than  the  DA  profile  and  there  are  significant  improvements  in 
behavior  of  the  intensities.  First,  the  streamwise  velocity  intensity  has  increased  and 
is  now  very  close  to  both  KMM  and  experiment.  Secondly,  the  artifact  produced  near 
the  wall  in  the  wall-normal  velocity  component  has  now  vanished  due  to  the  Green 
function  correction  which  insures  enforcement  of  continuity  at  the  wall.  This  implies 
that  wall  normal  derivative  of  u  >  must  vanish  at  the  wall.  There  is  also  a  marginal 
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increase'  in  the  wall-normal  velocity  intensity  and  a  somewhat  larger  decrease  in  the 
spanwise  velocity  intensity. 

In  figure  5  we  compare  the  Reynolds  stress  computed  from  flit'  DA  data  with  the 
KMM  data  and  the  experimental  data  of  Eckelmannf  1974)  who  made  measurements  at 
R*  values  of  142  and  20S.  The  DA  results  (recall  that  these  calculations  were  performed 
at  R*  =  12o)  are  in  reasonable  agreement  with  the  Eckelmann  (R*  =  142)  rt'sults  but 
are  considerably  lower  than  the  KMM  results  (R*  =  ISO)  and  the  higher  Reynolds 
number  experimental  values.  In  Figure  G  the  DAG  results  show  somewhat  better 
agreement  with  the  lower  Reynolds  number  experimental  results.  This  comparison 
suggests  that  for  the  relatively  weak  turbulence  represented  in  the  DA  and  the  KMM 
data  that  some  of  the  discrepancy  between  these  calculations  may  be  due  to  Reynolds 
number  differences. 

From  this  point  on  (Figues  7-1G)  we  will  make  comparisons  only  to  the  DAG 
datasets.  In  figure  7  the  correlation  coefficient  (Reynolds  stress  normalized  by  the 
relevant  root  mean  square  amplitudes)  shows  the  same  sharp  peak  near  the  wall  as  in 
the  KMM  data  although  its  value  is  somewhat  higher.  This  is  clearly  a  reflection  of 
the  lower  turblence  intensities  in  our  calculation.  We  note  however  that  the  peak  near 
the  wall  is  absent  in  the  experiments  of  Sabot  and  Comte-Belot  (197G). 

Tht'  skewness  for  each  velocity  component  is  plotted  in  Figures  S,9  and  10.  The 
data  generally  agree  with  the  KMM  results  anti  experimental  data  of  Barlow  and  John¬ 
ston  ( 1985)  and  Kreplin  and  Eckelmann  (1979).  The  notable  exception  is  flit*  normal 
component  of  velocity  where  KMM  and  DAG  agree  but  differ  significantly  with  the 
experimental  data.  The  flatness  for  each  velocity  component  is  plotted  in  Figures  11. 
12  find  13.  The  agreement  between  t lit'  DAG  and  the  experimental  data  is  again  juite 
good. 

The  root  mean  square  vorticity  for  DAG  is  plotted  Figure  14  and  shows  excellent 
agreement  with  the  KMM  data  and  the  experiments  of  Kastrinakis  anti  Eckelmann 
!  19S3 I  he  DAG  spanwise  vorticity  shows  precisely  the  same  minimum  near  the  wall 
at  about  ./•'  --  7)  as  in  the  KMM  result.  KMM  postulate  the  existence  of  wall  layer 
vortical  structures  to  explain  this  unexpected  result.  The  flatness  and  skewness  for  the' 
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fluctuating  component  of  vorticity  are  plotted  in  Figures  15  and  16  for  which  there  are 
no  comparisons  available. 

Generally  it  can  be  concluded  from  these  comparisons  that  the  mean  statistics  of 
the  DAG  dataset,  compare  well  with  the  highly  resolved  data  of  KMM.  Our  calculations 
confirm  the  minimum  in  the  spanwise  vorticity  as  seen  in  the  KMM  work.  These  cal¬ 
culations  also  confirm,  in  agreement  with  KMM,  spanwise  and  wall-normal  turbulence 
intensities  that  are  lower  than  experimental  values. 

It  is  clear  that  the  use  of  the  Green  function  method  improves  the  one-point 
statistics  computed  with  the  unmodified  Orszag-Kells  method  (i.e.  the  DA  dataset). 
However,  an  important  issue  has  been  raised  by  the  results  presented  above.  This 
issue  stems  from  the  observation  that  simulations  like  those  discussed  here  produce 
one-point  statistics  which  are  remarkably  close  to  those  produced  by  the  fully  resolved 
simulations  of  KMM.  It  certainly  cannot  be  claimed  that  if  one  achieves  good  one- 
point  statistics  that  fundemental  turbulent  dynamics  are  also  being  captured.  Indeed, 
conditional  sampling  performed  on  the  current  data  (Section  3.3)  indicates  a  bursting 
rate  that  is  lower  than  experiment.  On  the  other  hand,  it  is  difficult  not  to  conclude 
that  the  underlying  large  scale  (or  coherent  structures),  which  are  being  adequately 
resolved  at  the  resolutions  presented  here,  must  be  dominating  the  global  behavior  of 
the  flow  so  that  the  modest  resolutions  employed  here  are  adequate  to  achieve  good 
overall  results.  Certainly  this  is  not  a  strong  argument  to  justify  the  use  of  marginal 
resolution  but  it  does  suggest  that  expansion  functions  determined  from,  for  example, 
a  Karhunen-Loeve  expansion  (Lumley,1970)  may  be  adequate  to  capture  the  essential 
dynamics. 
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3.2  Two  point  statistics 


Using  the  DA  data  set,  we  have  computed  the  complete  two  point  correlation 
function  for  the  flow  in  the  following  manner.  For  any  given  realization  of  the  flow, 
u(x.f),  we  compute  its  Fourier  coefficients,  by: 


A  1  —  1  3  —  1 

2)=  ^  y;  u(/,  k,  X2  )e 
1=0  k= 0 


2  t  n  I  i 

*1 


2  it  tn  ft  1 

e  n3 


where 


i  =  x/^T, 
n  =  0, ....—  h  1, 


and 


(3.5) 


In  these  expressions  we  have  implicitly  defined  collocation  points  (x'i)i  and  (x3)k,  and 
wavenumbers  (ki)n  and  (£3),,,  as  follows: 


( x  1 )/  =  /  x  Xi/A'i, 


(,r3)A.  =  fc  x  L3/iV3, 


Z7T 

(A'i  )„  =  n  X  — , 
■01 

2?r 

( ^’3  )  m  =  X  — , 

■03 


(3.6) 


in  which  Li  and  Z,3  are  the  domain  lengths  in  the  streamwise  and  spanwise  directions 
respectively.  In  (3.5)  time  serves  merely  as  a  parameter  which  may  be  used  to  identify 
any  given  flow  realization  and  is  therefore  suppressed.  We  now  follow  Sirovich(19S7) 
who  suggests  that  the  geometric  symmetries  of  a  flow  should  be  exploited  to  increase 
the  effective  number  of  realizations  that  are  used  in  the  averaging  process.  For  the  case 
of  channel  flow,  the  equations  of  motion  are  invariant  with  respect  to  wall  normal  re¬ 
flection.  spanwise  reflection,  and  rotation  about  the  streamwise  axis.  These  symmetries 
yield  additional  spectral  representations  of  the  flow  which  can  be  written: 
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=  {$i0(n,m,-j2),-$2°(n,in,-X2),$30(n,in,-i2)} 

3r  =  {$i°(n,-m,x2),$20(n,-m,x2),-<I>30(n, -m,:r2)} 

#3  =  {$i°(n,  in,  -x2),  -$2°(n,  -m,  -x2),  -$3°(«,  -m,  -x2)}. 

(3.7) 


In  these  expressions,  the  superscripts  1,  2,  and  3  refer  to  vertical  reflection,  span- 
wise  reflection,  and  streamwise  rotation  symmetries  respectively  and  the  zero  super¬ 
script  refers  to  the  unaltered  (spectral)  representation  of  the  flow.  The  subscripts 
continue  to  represent  the  coordinate  directions  as  previously  defined.  The  (complex) 
spectrum  'I' Qa  can  now  be  computed  from: 

l  3  _ 

Va0{n,m,x2,x2')  =  (-J '2$Q,,(n,m,x2)$ljp(n,m,x2 '))  (3.8) 

p= o 

where, 


a,  ft  =  1,2,3. 

The  overbar  designates  complex  conjugate,  and  the  brackets  represents  an  average  over 
all  30  realizations  in  the  DA  dataset.  We  note  that  in  addition  to  the  increase  in  the 
effective  size  of  the  data  base  by  a  factor  of  four,  the  symmetries  have  the  additional 
advantage  of  reducing  the  size  of  'kQjj  also  by  a  factor  of  four.  A  factor  of  two  comes 
from  the  spanwise  symmetry  which  insures  that  is  unique  for  n  >  0  and  m  >  0. 
A  second  factor  of  two  comes  from  the  vertical  reflection  symmetry  which  gives  unique 
values  of  the  spectrum  for  —  1  <  x2  <  0  and  —  1  <  x2  <  1.  The  two  point  correlation 
function  ,  RQ/j,  can  now  be  obtained  from: 


2 


"a_i 

2 


R, 


,jU,  k,x2,x2)  =  K  ^  ^  ^Qt){n,m,x2,x2')e 


2 *  ii  l  i  2ym<fi 
N,  f  ,V3 
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n  —  (l 


2 
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resolution  remaining  unchanged)  and  also  by  running  at  higher  absolute  grid  resolution. 
We  also  note,  in  Figure  ISb,  a  clear  example  of  the  errors  induced  by  the  Orszag-Kells 
algorithm  near  the  wall.  The  curve  labeled  a ,  for  r2*  =  0.602  appears  to  be  completely 
unrelated  to  the  other  curves  in  Figure  18b.  For  this  distance  from  the  wall  the  vertical 
component  of  velocity  undergoes  a  much  too  rapid  loss  of  streamwise  correlation.  This 
peculiar  behavior  must  be  connected  to  the  boundary  condition  inconsistency  problem 
described  in  Section  2.  It  is  interesting  to  note,  however,  that  this  error  shows  up  in 
the  correlation  function  only  for  the  vertical  velocity  component  and  then  only  very 
near  the  wall.  This  does  not  imply  however,  that,  other  errors  are  produced  which  are 
not  localized  to  the  wall  region.  Indeed,  as  described  in  Section  3.1,  we  found  some  loss 
of  velocity  fluctuation  intensity  in  the  buffer  region.  Finally,  we  note  the  interesting 
behavior  in  the  streamwise  correlation  curves  for  the  spanwise  velocity  shown  in  Figure 
18c.  The  clear  separation  of  the  curves  below  ,r2*  =  14.76  and  those  above  ,r2*  =  28.3 
is  evident  which  suggests  that  the  spanwise  velocity  structure  is  particularly  sensitive 
to  the  differences  between  inner  and  outer  layer  physics. 

In  Figure  19  we  present  the  wall  normal  dependence  of  the  correlation  function 
by  plotting  i?Q/j(x2,  ,r2').  That  is  we  plot  Ralj  given  by  (3.9)  with  l  =  k  —  0.  The 
correlations  for  the  streamwise  velocity  (Figure  19a)  are  positive  except  for  a  weak 
negative  correlation  for  45.7  <  .r2*  <  77.2  which  appear  at  values  of  r2'  across  the 
centerline  of  the  channel.  Such  weakly  negative  correlations  also  appear  in  isotropic 
turbulence  (see  Hinze,  1975).  The  wall  normal  velocity  correlations  shown  in  Figure 
19b  are  positive  everywhere  with  the  exception  of  the  wildly  oscillating  positive  and 
negative  behavior  for  x2*  =  0.602.  Again,  this  behavior  is  attributed  to  the  operator 
splitting  errors  discussed  previously.  In  Figure  19d  we  present  the  results  for  i?]2  and 
note  that  this  function  does  not  peak  when  x2'  =  ,r2  and  indeed  there  is  no  theoretical 
reason  for  the  peak  to  occur  there.  We  note,  for  example,  that  close  to  the  wall,  say 
for  0  <  ,r2'*  <  25.  the  peaks  occur  clearly  at  values  of  ,r2'*  for  which  x2'*  >  x2*. 

The  spanwise  velocity  correlation  function  shown  in  Figure  19c  shows  the  most 
interesting  behavior.  We  first  note  that  for  the  curves  0  <  ,r2*  <  21.07,  the  correlation 
is  always  positive  as  we  move  into  the  wall.  The  curve  for  ,r2*  =  28.37  shows  the 
first  negative  correlation  values  near  the  wall.  All  curves  show  negative  correlation  for 
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some  values  of  x2'.  when  x2'  >  ■  One  possible  interpretation  of  these  correlations 

is  that  near  the  wall  there  is  a  streaxnwise  oriented  vortex  at.  x>*  ~  20.  That  is,  we 
interpret  the  curve  for  x2 *  =  21.07  as  the  center  of  the  wall  vortex  since  its  correla¬ 
tion  coefficient  is  almost  zero  hilt  still  slightly  positive  near  the  wall.  This  location 
of  the  wall  vortex  center  is  consistent  with  the  interpretation  of  the  vorticity  intensity 
profiles.  The  second  minima  (i.e.,  the  maximum  negative  values  which  exist  for  each 
curve  for  ,*•>'  >  x2)  moves  away  from  the  wall  roughly  in  proportion  to  x2.  That  is. 
the  scale  of  motion  seen  by  observers  farther  from  the  wall  is  proportional  to  their 
distance  from  the  wall.  We  remark,  however,  that  any  attempt  to  interpret  these,  or 
any  other  long-time  correlation  functions  in  terms  of  some  proposed  flow  structures 
without  the  addition  of  flow  visualization  information  is  highly  speculative.  That  is. 
there  could  be  many  possible  flow  structures  which  could  be  used  to  explain  any  set 
of  correlation  functions.  The  striking  features  of  the  R33  function,  however,  do  ap¬ 
pear  to  be  in  agreement  with  the  currently  accepted  viewpoint  that  streamwise  vortex 
structures  exist  near  the  wall.  We  also  remark  that  perhaps  the  most  objective  way 
of  reconstructing  a  given  three  dimensional  flow  from  the  complete  three  dimensional 
correlation  function  is  by  means  of  more  sophisticated  statistical  approaches  such  as 
the  orthogonal  decomposition  proposed  by  Lumleyf  1070). 

Finally,  in  Figure  20  we  compare  our  calculations  with  the  measurements  of  C’omte- 
Bellot(  1965).  We  note  first  that  these  measurements  were  at  a  Reynolds  number  of 
30.000.  a  factor  of  12  greater  than  in  our  calculations.  This  accounts  for  the  consid¬ 
erably  more  rapid  decay  in  the  measured  correlations.  We  note  in  particular  that  the 
measurements  show  small  negative  correlation  values  for  J?n  for  the  middle  curve  of 
Figure  20a  and  the  two  (negative)  minima  in  R 33.  These  characteristics  are  both  in 
excellent  agreement  with  our  calculations. 
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3.3  Conditional  Sampling  Results 


3.3.1  Purpose 

Conditional  sampling  of  a  turbulent  channel  flow  is  a  heuristic  search  for  events 
or  conditions  which  appear  to  have  some  relevance  to  the  kinematics  and  dynamics  of 
the  flow.  The  search  is  usually  based  on  some  prc-conceived  idea  or  hypothesis  which 
is  used  as  a  condition.  When  the  flow  is  conditionally  sampled,  the  statistics  about 
the  occurences  of  the  events  can  be  determined.  Some  statistics  of  interest  may  be  the 
frequency  of  occurence  of  the  detected  events,  the  ensemble  average  of  the  events,  or  a 
histogram  of  some  property  associated  with  the  detected  event. 

Conditional  sampling  has  been  used  by  many  investigators  to  detect  turbulent 
bursts.  A  turbulent  burst  is  an  ill-defined  event  or  sequence  of  events  which  is  re¬ 
sponsible  for  the  generation  of  a  disproportionate  amount  of  the  Reynolds  stress  and 
turbulent  kinetic  energy.  In  a  single  turbulent  burst,  there  may  be  several  ejections 
of  low  speed  fluid  from  the  wall  region  into  the  core.  In  addition,  associated  with  the 
turbulent  burst  is  a  shear  layer  which  is  usually  detected  as  a  sweep  (1*2  >  0,  <  0). 

The  ensemble  average  of  the  conditionally  sampled  events  is  considered  to  be  a 
ropresentive  event  of  the  flow  physics  for  a  short  period  of  time  (say  10f*).  The  rapid 
fluctuations  seen  in  turbulent  flow  limit  the  significance  of  the  ensemble  after  a  short 
time. 

3.3.2  The  One  And  Two  Point  Detections 

There  were  two  objectives  to  the  one  and  two  point  detection  tests.  The  first 
objective  was  to  determine  if  similar  timing  statistics  and  ensemble  averages  could  be 
obtained  from  the  numerical  data  set.  Two  standard  detection  schemes,  the  Variable 
Interval  Time  Averaging  scheme  (VITA)  and  the  Second  Quadrant  Detection  scheme 
( Q> )  were  used.  They  were  applied  without  taking  advantage  of  the  spatial  character 
of  the  data:  that  is  they  were  used  as  an  experimentalist  would  have  applied  them. 
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The  second  objective  was  to  determine  the  proximity  of  two  events  in  space  and  time. 
For  example,  a  particular  test  was  to  determine  how  close  two  related  Q 2  events  are  in 
time. 

The  VITA  detection  scheme  is  based  on  the  observation  that  large  fluctuations  in 
the  streamwise  velocity  occur  during  a  typical  burst.  These  large  fluctuations  are  caused 
by  a  shear  layer  passing  by  the  detector.  The  VITA  technique  measures  the  strength 
and  thickness  of  the  shear  layer  by  using  a  threshold  parameter,  k,  and  an  averaging 
time  interval,  Ta.  When  the  short  time  variance,  properly  normalized,  exceeds  the 
threshold,  k,  a  burst  is  detected.  A11  additional  condition  of  the  sign  of  frequently 

used.  Unless  noted  in  this  work  only  positive  slope  events  >  0)  are  retained. 

The  Reynolds  stress  can  be  sorted  into  the  four  quadrants  in  the  u\  —  u'2  plane. 
(See  table  3.2)  The  second  quadrant  Reynolds  stress  algorithm,  or  the  Q2  detection 
algorithm  is  based  on  the  observation  that  to  be  important,  turbulent  bursts  must 
develop  Reynolds  stresses,  and  second  quadrant  Rejmolds  stresses  would  be  the  defin¬ 
ing  property  of  turbulent  ejections.  In  this  detection  scheme,  an  event  is  considered 
significant  if  the  magnitude  of  second  quadrant  Reynolds  stress  rises  above  a  predeter¬ 
mined  level.  Both  of  these  algorithms,  along  with  others,  are  described  in  Luchik  and 
Tiedermann  (1987). 

Table  3.2 

Reynolds  Stress  Quadrant 


Quadrant 

w'i 

u'2 

Physical  Process 

Qi 

>  0 

>  0 

Interaction 

Q2 

<  0 

>  0 

Ejection 

Q 3 

>  0 

<  0 

Interaction 

Q-x 

>  0 

<  0 

Sweep 
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Table  3.3 

Intra-Burst  Ejection  Time  Scale 


UJevel  Q2  threshold. h  Intra-Burst  ejection  timescale 


0.0 

0.5 

12.3 

0.0 

1.25 

15.2 

0.0 

2.0 

19.9 

-0.5 

0.5 

11.2 

-0.5 

1.25 

13.0 

-0.5 

2.0 

1G.0 

-1.0 

0.5 

9.9 

-1.0 

1.25 

11.0 

-1.0 

2.0 

11.9 

-1.5 

0,5 

9.5 

-1.5 

1.25 

10.4 

-1.5 

2.0 

10.6 

The  timescale  discussed  above  relates  two  ejection  events  passing  a  single  detector. 
Of  related  interest  is  the  probability  of  a  second  event  being  detected  after  an  initial 
detection.  For  example,  how  likely  is  the  detection  of  a  Q2  event  at  a  location  .Ti  and 
time  T2  given  that  a  detection  has  occured  at  a  different  location  x\  and  time  7j .  We 
may  also  ask  how  likely  is  a  VITA  event  to  follow  or  lead  a  Q 2  event  by  some  time 
span.  With  this  approach  the  relationship  between  differing  detection  schemes  can  be 
studied. 

The  simulated  flow  was  conditionally  sampling  at  32  points  in  the  spanwise  direc¬ 
tion  ( A:r*  =  20)  at  two  streamwise  locations.  Then  the  following  probability  variable 
was  developed.  It  is  not  a  normalized  probability,  but  rather  a  probability  relative 
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to  a  random  occurence.  Before  defining  the  probability  variable,  let  us  introduce  the 
following  notation:  Ti(i ,  n)  and  T2(i.n)  is  the  detection  time  of  an  event  i  at  spanwise 
location  n  and  streamwise  location  a* i ( 1 )  and  J'i(2),  respectively,  where  x  i  ( 1 )  is  the 
primary  location  and  Xi(2)  is  the  secondary  location.  Furthermore,  introduce  a  time 
shift  T<  =  T>  —  T\ .  The  parameter  we  want  to  define  is  the  probability  of  detecting  an 
event  i  at  the  location  T>(i,n)  given  a  detection  j  at  Ti(j,  n).  This  parameter  is  there¬ 
fore  a  function  of  the  two  detection  schemes  applied,  the  shift  in  time,  T,,  and  the  shift 
in  the  streamwise  direction  X3  =  .r i ( 2 )  -  .r,(  1 ).  We  therefore  define  the  probability,  p. 
as  follows: 


p{  -V, .  J, )  = 


number  of  events  at  T2(i,n)  given  a  detection  at  Ti(?\n)}, 


where  is  the  total  number  of  events  detected  at  streamwise  location  .r  j  ( 1 ).  For  large 
time  separations,  the  events  at  the  two  streamwise  locations  become  uncorrelated,  and 
the  probability  becomes  the  frequency  of  the  detection  of  the  T>  events.  Therefore,  if 
the  probability  above  is  renormalized  by  the  mean  detection  period  of  the  T2  events,  the 
resulting  probability  is  a  measure  of  the  likelihood  that  a  detection  will  occur  relative 
to  the  random  detection. 

Some  of  the  main  results  are  as  follows: 

(1)  Positive  slope  VITA  events  are  less  dispersive  as  they  convect.  When  the  primary 
and  secondary  detectors  use  VITA  and  A”*  =  240,  two  peaks  in  the  relative  probability 
can  be  seen.  The  first  peak,  centered  in  time  about  the  convective  shift  20f*.  is  a 
result  of  those  events  detected  at  the  primary  location.  A  second  peak,  centered  about 
a  convective  shift  of  — 40f*  is  a  result  of  those  events  which  will  be  detected  at  the 
primary  location. 

(2)  For  A*  —  240,  the  Q2  events  show  a  similar  two  peaked  behavior,  but  the  peaks 
have  spread  wider  than  for  the  VITA  events.  This  implies  that  the  convection  speed 
for  the  Q>  events  is  more  varied  than  for  the  VITA  events. 
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Table  3.4 

Program  Flow  Chart  for  Detection  and  Indexing  of  Streaks 

Input /Output 

Routine 

Process 

Velocity  Fields 

— >  Lss.ftn 

1 

< —  Resort.  Ft  n 

1 

< —  Timing.  Ftn 

Finds  All  Streaks 

Sorted  Streaks 

Sorts  Streaks  In  Time 

Sorted  and  Indexed 

Streaks 

Develops  Statistics 

Deletes  Short  Streaks 

Typical  low  speed  streaks  can  be  seen  in  Figures  23  and  24.  Figure  23  is  a  plan 
view  of  the  plane  y*  —  15.  Figure  24  is  a  view  looking  upstream  at  the  point  x*  =  0. 
Both  figures  are  of  the  strearnwise  velocity.  The  shaded  regions  indicate  where  the 
fluctuating  velocity  is  less  than  — u'i.  The  markers  indicate  where  the  streaks  have 
been  detected  and  followed.  There  are  ten  detected  streaks  at  this  point  in  time. 
Three  in  the  vicinity  of  x*  =  400  and  —  300  are  weakly  defined. 

The  spacing  of  the  detected  streaks  is  about  100  to  130  viscous  lengths  apart, 
which  is  consistent  with  experiments.  Depending  on  the  parameters  used  in  the  Lss.ftn 
program,  the  average  streak  spacing  was  between  105  and  115  viscous  lengths. 

We  also  find  (Fig.  24)  that,  the  locations  detected  occur  under  or  near  ‘tunnel1 
shaped  shear  surfaces.  Consider  the  detections  at  .r*  —  90  and  190,  which  were  typical 
of  smooth,  steady,  non-meandering  low  speed  streaks.  The  locations  =  340  and 
530  come  after  the  passage  of  a  streak  and  are  weakly  defined  since  the  streaks  may 
have  lifted  away  from  the  wall.  Note  that  inflection  velocity  profiles  can  be  seen  near 
./■*  =  70.  350  and  520.  all  of  which  are  near  low  speed  streaks. 
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Figure  25  shows  the  evolution  of  the  normal  velocity  along  a  low  speed  streak  in 
the  streamwise  direction  (the  horizontal  axis)  as  a  function  of  time  (the  vertical  axis). 
The  convection  speed  of  the  velocity  perturbations  remains  about  10a*  even  though 
the  streamwise  velocity  in  the  streak  is  about  S  —  9a*.  A  typical  pair  of  Q2  events  is 
located  at  t*  —  50  and  =  200.  At  t*  =  20  and  x*  =  350  a  perturbation  in  the  normal 
velocity  is  developing  in  front  of  an  existing  one.  This  is  the  development  of  a  new  Q> 
event.  The  existing  and  the  new  perturbation  may  be  related  to  the  pair  at  50,200. 

Unfortunately,  the  scheme  which  finds  the  initial  streaks  is  not  very  robust.  If  a 
streak  develops  a  kink  or  dips  below  the  plane  in  which  the  searching  is  done,  it  may  be 
evaluated  as  two  streaks  and  the  continuity  in  time  is  lost.  Experimentally,  the  effects 
of  low  speed  streaks  can  be  seen  for  a  few  thousand  viscous  time  units,  but  when  using 
the  algorithms  developed  on  the  data  available,  less  that  3%  of  the  detected  streaks 
lasted  more  than  200f*.  Better  results  are  expected  when  the  algorithm  is  applied 
closer  to  the  wall,  where  the  flow  is  more  stable. 
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4.  Conclusions 

Wo  have  investigated,  in  some  detail,  t lie  statistical  properties  of  turbulent  flow  in 
a  channel  using  data  generated  from  moderate  resolution  direct  numerical  simulations. 
Both  long-time  averages  and  conditional  averages  were  computed.  We  have  concluded 
that  the  use  of  the  Green  function  method  results  in  turbulence  statistics  which  are 
significantly  closer  to  experiments  than  those  produced  with  the  Orszag-Kells  algo¬ 
rithm  t  Section  3.1).  However,  the  velocity  fluctuation  intensities  are  still  somewhat 
lower  than  those  produced  by  the  fully  resolved  K.MM  calculation.  This  is  particularly 
true  for  the  wall-normal  and  spanwise  velocities.  However,  the  vorticity  profiles  are 
in  excellent  agreement  with  the  KMM  results.  For  the  two-point  statistics  presented 
in  Section  3.2  we  find  that  the  spanwise  correlation  lengths  are  in  excellent  agreement 
with  the  KMM  results  when  Reynolds  number  differences  are  taken  into  account.  The 
streamwise  correlation  lengths,  particularly  those  of  the  streamwise  velocity  compo¬ 
nent.  are  somewhat  larger  than  the  KMM  result.  We  speculate  that  Reynolds  number 
differences,  box  length  differences,  and  resolution  may  all  be  contributing  to  these 
discrepancies.  The  wall-normal  dependence  of  the  spanwise  velocity  correlation  func¬ 
tion.  i?  j  ?  (  r  >•  3*2 1  )■  shows  behavior  which  is  suggestive  of  near  wall  steamwise  vortical 
structures.  The  wall-normal  dependence  of  each  correlation  function  shows  excellent 
qualitative  agreement  with  the  measurements  of  Comte-Bellot. 

The  three  conclusions  to  the  conditional  sampling  work  are:  (1)  the  low  speed 
streak  tracing  is  able  to  follow  in  space  and  time  the  low  speed  streaks,  but  the  method 
needs  to  be  applied  closer  to  the  wall:  (2)  to  determine  temporal  statistics,  a  dataset 
with  a  longer  domain  is  needed,  and:  (3)  the  results  from  the  two  point  conditional 
sampling  scheme  indicate  that  there  is  a  very  definite  correlation  between  VITA  and 
Q>  events  and  that  the  model  proposed  in  Leighton  (19SG)  is  at  least  partially  correct. 
Further  developments  in  conditional  sampling  will  require  an  improved  dataset  and  will 
be  finished  in  the  near  future. 

We  have  found  that,  aliasing  errors  have  significant  effects  on  the  final  steady-state 
solutions  to  the  Xavier-Stokes  equations  for  channel  flow  (Appendix  A).  When  the 
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calculation  has  converged  to  the  predicted  steady-state  value  of  the  wall  shear  stress  and 
aliasing  errors  are  not  eliminated,  we  find  that  both  the  mean  velocity  and  the  centerline 
velocity  converge  to  values  which  are  significantly  lower  than  experimental  values.  In 
addition,  the  velocity  fluctuation  intensities,  particularly  that  of  the  streamwise  velocity 
component,  are  too  low  in  amplitude  and  peak  at  distances  too  far  from  the  wall.  We 
conclude  that  aliasing  errors  produce  solutions  which  appear  stable  but  which  are  very 
much  like  the  flow  of  underdeveloped  turbulence.  Aliasing  errors  appear  to  have  a 
damping  effect  on  the  turbulence.  The  strength  of  these  effects  may  well  be  due  to  the 
low  resolution  of  these  calculations,  but  what  appears  to  us  to  be  most  interesting  is  that 
when  these  errors  are  removed,  the  solution  goes  quickly  back  to  that  of  experimentally 
observed  turbulence. 
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Appendix 

A.  Effects  of  Aliasing  Errors  on  the  Behavior  of  Channel  Flow  Turbulence  Calculations 

1.  Introduction 

It  has  been  known  for  some  time  that  aliasing  errors  can  have  unwanted  effects 

/• 

on  spectral  calculations.  In  fact,  in  low  resolution  calculations,  aliasing  errors  can  lead 
to  catastrophic  numerical  instabilities.  Before  discussing  the  details  of  the  effects  of 
aliasing  in  turbulent  channel  flow  calculations,  however,  we  first  give  a  brief  review  of 
the  mathematical  origin  of  aliasing  errors.  We  first  begin  with  a  simple  example  of 
aliasing  effects  arising  in  the  representation  of  a  continuous  one-dimensional  function. 
We  assume  that  we  are  given  a  continuous  function  which  is  represented  by  an  infinite 
(convergent)  Fourier  series.  That  is,  u(.r)  is  given  by: 

OO 

u(x)  =  ^  Uke'kx  (Al) 

k=—oo 


where 


dx. 


Now.  we  suppose  that  we  sample  the  continuous  function  at  the  discrete  points,  xj  = 
=  0....N  —  1.  It  follows  that: 


and. 


A'/ 2  - 

"(*•>)=  E 

k—  —  S/2 


uL.e 


ikxj 


N-l 


>=o 


(A2) 


(A3) 
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That  is,  we  clearly  distinguish  the  continuous  Fourier  modes  from  the  discrete  modes. 
To  find  the  relationship  between  these  two  sets  of  modes  we  write: 


A’  — 1  oo 


£  a*'*"'-*"- 


(A4) 


j= 0  k=  —  oo 


or 


oo  1  jV  —  1 


p  ^  "'EY  . 

k=  —  oo  j = 0 


where 


J_  ei(k-p)xj  _  f  1,  if  k  -  p  =  m.V  ?7?  =  0,  ±1. 

Ar  1  0  otherwise. 


(A5) 


Therefore,  we  find  the  following  relationship  between  the  discrete  and  continuous 
modes: 


OO 

u*  =  iip  +  Up+mN-  (AQ) 

m  r=  —  oo 
myt  0 

The  second  term  on  the  right  represents  the  aliasing  errors  which  contaminate  the 
discrete  Fourier  representation  of  u( x).  We  note  that  if  there  is  no  energy  in  u(x)  beyond 
mode  N  then  these  errors  vanish  and  the  discrete  coefficients  equal  the  continuous 
coefficients. 

In  Navier-Stokes  calculations,  aliasing  errors  arise  from  the  non-linear  terms  in  the 
equations  of  motion.  For  example,  the  term  u  x  u>  gives  rise  to  errors  when  the  rota¬ 
tional  form  of  the  equations  are  used.  In  pseudo-spectral  calculations  we  perform  the 
nonlinear  product  of  velocity  and  vorticity  in  physical  space  as  opposed  to  computing 
the  convolution  sums  which  would  arise  if  a  purely  Galerkin  procedure  were  used.  This 
procedure  of  performing  the  non-linear  calculation  in  physical  space  has  been  shown 
by  Orszag(  19711))  to  be  considerably  more  efficient  then  the  Galerkin  approach  but 
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then  aliasing  errors  may  contaminate  the  calculation.  A  simple  example  will  clarify 
the  mathematical  origin  of  aliasing  due  to  nonlinearities.  We  are  given  the  discrete 
representation  of  two  functions  u(x)  and  v(x )  as  follows: 


and 


where 


N/ 2-1 

u(xj)  =  utc'*Xi» 

k=-N/2 


N/2-1 


;(.rJ)=  5Z  ul 


P=-N/2 


0,1, . N  -  1. 


(-47) 


(A8) 


The  product  of  these  functions  is  then  given  by: 


k  p 


v/ei(k+p)x>, 


where 


Zj  =  u(Xj)v(Xj) 


and  the  Fourier  coefficients  of  zj  are: 


(A9) 


z 


* 


N-l 

__  ^  e(k+p-q)zj 
^  J= 0 


M10) 


Now  we  use,  again,  the  identity: 
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(All) 


TV— 1 

—  ei(k+p~q)x>  =  (  1’  k+p-q=mN  m  =  0,±1,...; 
N  \  0  otherwise. 

3=0 


(A  10)  becomes: 


OO 

zg  =  E  Uk*<  +  E  E  Uk*vP*' 

(k,p)  m  =  —  oo  (fc.p) 

fc  +  p  =  g  m^O  A;  -}-  p  —  q  =  m  y 


However,  we  note  the  following  inequalities; 


(A12) 


and 


|p U  HU  I  < 


N 
o  ’ 


(A13) 


max(p  4-  A;  —  q)  = 

Since  ^  <  2iV  it  therefore  follows  that  |  m  |<  1.  The  final  result  becomes: 

<=Ew**v  +  E  u**v+  E  «**v-  u^) 

(*,p)  (*.P)  (*.P) 

k+P  —  <l  k+p  =  q  +  N  k  +  p  =  q-S 

The  last  two  terms  on  the  right  hand  side  of  (A14)  are  termed  aliasing  errors.  We 
observe  that  if  all  modes  y  >  (k,p)  >  andyp  >  ( k.p )  >  -y  are  removed  from 
the  spectral  representation  of  u  and  v  and  if  their  product  is  taken  in  physical  space, 

then  there  can  be  no  aliasing  errors  in  the  region  y  >  q  >  y,  andyy  >  q  >  yp. 

This  fact  was  first  observed  by  Orszag  (1971a).  We  note  that  this  so  called  ‘two- 
thirds  rule’  works  because,  for  example,  when  k  +  q  >  =~-  aliasing  errors  are  dumped 
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into  the  region  >  q  >  Since  this  region  of  the  spectrum  is  nulled  after  the 

non-linear  product  is  taken,  the  calculation  is  alias-free.  We  note  that  to  insure  an 
alias-free  calculation  one-third  of  the  modes  is  the  minimum  number  that  need  to  be 
removed  and  is  therefore  an  optimal  choice.  In  a  typical  three  dimensional  calculation 
the  dealiasing  is  traditionally  performed  only  in  the  periodic  plane  and  therefore  in  our 
calculations  dealiasing  is  performed  using  the  ‘“2/3  rule”  in  the  x  i  —  X3  plane.  There 
is  some  evidence,  however,  that  in  the  non-homogeneous  direction,  (X'> )  that  aliasing 
effects  may  be  important  (Zang  and  Krist(1987)). 

2.  Results  from  low  Resolution  Aliased  Channel  Flow  calculations 

A  relatively  long  run  was  made  using  the  Orszag- Kells  algorithm  with  the  boundary 
condition,  112  =  0  on  X2  —  ±1,  applied  in  the  pressure  step  in  which  aliasing  errors 
were  not  removed.  We  used  33  x  64  x  16  grid  points  in  the  X2,3'3,and  xi  directions 
respectively.  The  calculation  was  started  with  laminar  initial  conditions  to  which  2D 
and  3D  finite  amplitude  Orr-Sommerfeld  modes  were  applied. 

Before  proceeding  with  a  detailed  description  of  the  results,  however,  we  review 
the  manner  in  which  the  flow  is  driven  to  a  steady-state.  This  is  discussed  in  Handler, 
Leighton,  and  Carroll(1987)  but  is  reiterated  here  for  the  sake  of  completeness.  First, 
we  note  that  in  a  physical  experiment  the  flow  is  driven  by  a  pressure  gradient.  In  the 
simulation,  however,  periodic  boundary  conditions  are  used  in  the  streamwise  direction 
so  that  a  differential  'pressure  forcing  cannot  be  prescribed.  Thus,  an  external  driving 
force,  f  .  is  required  in  the  equations  of  motion.  We  note,  of  course,  that  f  is  a 
completely  arbitrary  function  of  space  and  time.  If  we  define  a  control  surface  to  be 
that  which  encloses  the  entire  computational  domain,  a  momentum  balance  in  the 
streamwise  direction  gives: 


0_ 

Of 


<  Ui(x,t)  >  d.X2 


-2  <  rw  >  2  <  fi(x,t)  >  L2 

plr 2  +  pU2 


(-415) 
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In  this  expression,  <>,  represents  averaging  over  a  horizontal  plane,  tw  is  the 
instantaneous  shear  sueaa  at  the  wall,  and  f\  is  the  streainwise  component  of  the  force 
per  unit  volume.  The  simple  form  of  this  expression  is  due  to  the  use  of  periodic 
boundary  conditions  which  insures  no  net  mass  flow  through  the  control  surface  and 
no  net  forcing  arising  from  normal  stresses.  It  follows  that  for  a  globally  steady  flow 
the  time  rate  of  change  of  momentum  within  the  control  volume  must  be  zero,  so  that: 


lim  <  tw  >=  L2  lim  <  fi('x.J)  >  .  (.416) 

t — -OO  1  —  0 o 

This  can  be  written  in  non-dimensional  form  as: 

lim  R*(t)  =  lim  <  —  >  ^ 2  1 12.  (-417) 

i  —  OO  t  —  OC  pU  ~ 

Thus,  if  the  arbitrary  force  reaches  a  limit  as  t  — ►  o o,  we  can  predict  the  final  wall 
Reynolds  number. 

Some  global  results  of  this  aliased  calculation  are  presented  in  Figures  Al(a,  b,  c). 
In  this  calculation,  the  force,  -kjjf ,  was  fixed  at  ^  and  R  —  4000,  so  that  we  can  predict 
from  (A17)  that  the  final  Reynolds  umber,  i?*,  will  be  s/AR  =  126.49.  We  observe  in 
Figure  Al(a)  that  R *  initially  increases  rapidly  as  in  the  calculation  of  Orszag  and 
Patera  (1983)  but  then,  after  the  flow  has  gone  through  transition  to  turbulence,  the 
wall  Reynolds  number  gradually  decays  to  the  predicted  value.  Some  oscillations  about 
the  predicted  value  of  126.49  occur  but  a  dynamic  equilibrium,  that  is  a  state  in  which 
the  driving  force  is  balanced  by  shear  stresses  at  the  wall,  has  been  achieved.  We  can 
now  use  the  correlations  introduced  by  Dcan(1978)  for  fully  developed  channel  flow  to 
determine  values  for  Rci  and  Rm,  the  Reynolds  number  based  on  centerline  velocitj', 
Uri-  and  that  based  on  the  average,  or  bulk,  velocity  Um .  We  perform  this  calculation 
by  using  the  value  R*  =  126.49.  The  relationships: 

1-77-2  =  0-073R,„-25  (.418) 

2  Pt  ni 
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and 


r 


=  u 


U  =: 


R* 

R 


can  be  used  to  predict  i?m,  and  the  relation 


(419) 


(A20) 


=  1.28i?„r00116  (A21) 

U  m 

is  used  to  predict  the  centerline  velocity  which  can  then  be  used  to  compute  i?cj.  The 
results  of  this  calculation  give  Rm  —  3697.5  and  Rci  =  2151.3.  From  Figures  Al(b) 
and  Al(c)  we  find  that,  at  the  end  of  the  calculation,  Rm  =  4470  and  Rci  =  2800. 
Thus,  we  find  that  both  the  mean  Reynolds  number  and  centerline  Reynolds  number 
are  significantly  higher  then  they  would  be  in  an  experiment  with  this  asymptotic  value 
of  the  wall  shear  stress.  The  effect  of  aliasing  has  been  to  give  a  larger  mass  flux  for 
the  given  pressure  gradient  (driving  force)  than  for  fully  developed  turbulence.  This 
gives  a  strong  indication  that  the  flow  is  not  a  fully  developed  turbulence.  Evidence 
for  this  is  given  in  Figures  3(a,  b)  in  which  we  have  performed  a  somewhat  different 
calculation.  Here  we  have  used  a  fully  converged  solution  which  we  obtained  from 
Steven  Lyons  at  the  Univ.  of  Illinois  from  an  alias-free  code  as  an  initial  condition 
to  our  (aliased)  code.  We  note  that  the  effect  of  aliasing  was  first  (Figure  A2(a))  to 
increase  the  centerline  velocity  and  cause  a  dramatic  loss  of  the  logarithmic  region 
of  the  mean  velocity  profile  that  was  clearly  present  in  the  initial  condition.  The 
final  velocity  profile  is  very  similar  to  that  obtained  by  Patel  and  Head  (1969)  in 
transitional  (underdeveloped)  turbulent  channel  flow.  We  also  observe  in  Figure  A2(b) 
that  the  intensity  of  the  streamwise  velocity  fluctuations  are  dramatically  reduced  and 
the  location  of  the  maximum  intensity  has  shifted  farther  from  the  wall. 
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We  also  report  on  an  independent  test  of  the  effects  of  aliasing  performed  by  Panda 
(19S7)  who  performed  low  resolution  channel  flow  calculations  in  which  the  dealiasim 
algorithm  was  turned  on  at  an  appropriate  time  in  the  calculation.  These  results  ai< 
shown  in  Figures  A3(a)  and  A3(b).  We  note  that  the  effect  of  dealiasing  was  to  can-' 
a  dramatic  decrease  in  the  centerline  velocity  which  is  consistent  with  our  results 
There  can  be  little  doubt  from  these  results  that  aliasing  errors,  at  least  for  these  low 
resolution  calculations,  has  a  damping  effect  on  the  turbulence.  We  observe  first  that 
the  effect  is  not  clearly  exhibited  unless  the  calculations  are  run  for  very  long  rime-, 
so  that  these  errors  have  a  chance  to  build  and  interfere'  in  some  way  with  the  prope-i 
solution  to  the  equations.  Secondly,  and  perhaps  most  importantly,  the  aliasing  error- 
caused  no  catastrophic  numerical  instabilities.  Instead,  the  effects  were  subtle1  in  tin' 
sense  that  earlier  in  the  calculations  (results  not  shown  here)  the1  mean  velocity  profile 
did  exhibit  a  logarithmic  layer  which  subsequently  disappeared.  In  cemclusion.  aliasing 
errors  affect,  in  a  significant  manner,  the  long-time  behavior  e>f  the1  turbulent  solutie>n 
to  the  equations  of  motion.  We  speculate  that  aliasing  errors  are  acting  very  niue-h 
like  an  external  forcing  which  has  both  the  proper  phase  and  amplitude-  to  cause  some- 
damping.  but  not  a  complete  destruction  of  the  turbulence. 
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B.  Calculation  of  Energies  from  Fourier  Spectra 

When  Fourier  spectra  are  computed,  as  in  (3.8),  it  is  useful  to  check  the  results  by 
computing  from  them  the  mean  square  values,  or  energies,  by  integrating  the  spectra 
over  all  wavenumbers.  This  result  can  then  be  compared  with  a  'direct'  calculation  of 
the  energy  by  squaring  and  averaging  the  physical  variable  of  interest  such  as  velocity 
or  pressure,  for  example,  in  its  physical  (space-time)  representation.  We  derive  here  a 
convenient  expression  for  a  real  function  of  two  (space)  variables  that  can  be  used  to 
make  such  a  check.  We  first  define  the  real  function  u(z,y)  and  its  Fourier  coefficients 
ctk.p  as  follows: 


S/2-1  S/2-1 


U  l,i 


E  E 


[B 1) 


A- -A'/ 2  p  —  —  \'/2 


and. 


S-\  S- 1 


ak,p  =  T,u^e~,kr‘e~,Pym-  (B2) 

i-o  »i=0 

The  physical  lengths  are  defined  by  xj  =  ^j-l  and  ym  =  jj-m.  The  mean  square  energy, 
i/2.  is  given  bv  : 

(B3) 


S  —  1  N  —  1 

1,2  =  wEE ui^2- 


/=:  0  (11  =  0 


It  follows  from  (Bl)  that: 

^  =  E  E  m  E  E  c“+l'>‘vw>.. 


A'-l  A'-l 


(B  4) 


k,p  k' ./>' 


1=0  m  =  0 


Use  of  the  identities: 


A  —  ] 

—  =  J  L  lf  l>  +  p'  =  n.\  u  =  U.±l 

A’  y  l  0  otherwise. 

di=0 


(T?5) 
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and 


J_  y'  ei(k+k')t,  _  r  1,  if  k  +  k'  —  nN  n  =  0.  ±1, ...; 

N  1 0  otherwise. 

1=0 


(B  6) 


and  the  use  of  the  reality  of  u.  that  is  auP  =  a-k-P<  lead;,  to  the  following  expression 
for  the  energy: 


9  2,2  .  2  ,  2  , 
u~  ~  G0,0  +  a-.\'/2,0  w  °0.-.V/2  "T  a-A’/2.-.V/2  ' 

(V/2  — 1  jV/2-1  .V/2-1 

2{  51  i°*t>i2+  51 

t'=l  p=-.V/2+l  P=1 

JV/2-1  iV/2  — 1 

+  5]  |«  — A72,/>|2  +  5Z  A'/2  1“  }  - 

p=l  k  = 1 


(57) 


As  discussed  in  Section  3.b,  it  is  sometimes  useful  to  impose  the  symmetry  condition 
cik.p  —  cu,_p .  If  this  is  done  the  energy  becomes: 

"" 2  2  2  o 

U2  =  G0,0  +  a-.V/2,0  +  a  0,-.V/2  +  «1a72,-A72  + 

■V/2-1  A72-1 

4  53  53  la^i2 

fc=l  p=l 

A72-1  .V/2-1 

+2{  53  i«a-,o|2  +  53  io°/> i2 ) 

k  =  1  /»=  1 

.V/2-1  .V/2-1 

+2{  5^  la-A72.pl2  +  53  k .v/2 12 } ■  (58) 

p=l  A=1 

If  we  define  the  Fourier  coefficients  as  follows: 

«o,o  =  2.4o,u  a-jV/2.0  =  2A_  v/2,0 
a0.-.V/2  =  2-40._.V/2  a-N/2.-.\'/2  =  2-4  —  A/2.  —  .V/2 
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<H-,o  —  n/2.4/;, o  cio.p  —  v^-4o,P 

«-,V/2  .p  =  '/2-4_V/2,P  Uk,-N/2  =  V/2^1fc,-.V/2 

and. 

®k.p  —  -4jt,p 

1  <  A-  <  JV/2  -  -  1 
1  <  p  <  A*/2  -  1. 

then  the  energy  can  be  conveniently  expressed  as  : 

JV/2  Ar/2 

“2  =  45Z  £  I-4*-',pI2- 

r=o  p=o 

( B 1 1 )  was  used  to  compute  the  mean  square  energies  from  the  spectra,  4'aQ-  The 
results  for  the  root  mean  square  values  for  the  streamwise,  spanwise,  and  wall-normal 
turbulence  intensities  are  shown  in  Figure  B1  and  for  the  Reynolds  stress  ,ufu 2,  in 
Figure  B2.  One  can  see  that  these  results  are  identical  to  the  results  obtained  by 
performing  the  ‘direct’  computation  described  in  section  3.1  and  shown  in  Figures  2 
and  5.  We  conclude  with  a  high  degree  of  confidence  that  our  calculaton  of  is 
correct. 


(B  9) 


(510) 


(511) 
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Mean  Velocity  Profiles  from  the  DA  dataset.  □,  Current  results.  — .  KMM,  — ,  Hussain  and 
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Fig.  2  —  Velocity  Fluctuation  Profiles  from  the  DA  dataset.  □,  uu  (solid  triangle),  u2,  •,  u3.  Current 
results;  — ,  ut,  - — ,  u2, - ,  m3,  KMM;  □,  w,,  A,  u2,  o,  u3,  Kreplin  and  Eckelmann  (1979). 
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Fig.  4  —  Fluctuating  Velocity  Profiles  from  the  DAG  dataset.  □,  u  f,  (solid  triangle),  u2,  •,  u 3,  Current 
results;  — ,  ut,  — ,  u2,  —  -  «3,  KMM;  □,  uh  A,  u2,  o,  u3,  Kreplin  and  Eckelmann  (1979). 
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Current  results;  —  KMM; 
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Skewness  for  u2 ■  A,  Current  results;  — ,  KMM;  +  ,  Barlow  and  Johnston  (1985);  o,  Kreplin  and 
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Fig.  10  —  Skewness  for  u3.  •,  Current  results;  — ,  KMM;  o,  Kreplin  and  Eckelmann  (1979). 


Fig.  11  —  Flatness  for  u,.  o,  Kreplin  and  Eckelmann  (1979).  □.  Current  results;  — .  KMM;  +  ,  Barlow  and 
Johnston  (1985);  O,  Kreplin  and  Eckelmann  (1979). 


Flatness  for  u:.  A.  Current  results;  —  KMM:  +.  Barlow  and  Johnston  (1985);  o.  Kreplin  and 
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Fig.  13  —  Flatness  tor  u  t.  •,  Current  results;  — ,  KMM;  o,  Kreplin  anti  F.t  kclmann  (ll)71)). 
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Fluctuating  Vorticity  Profiles,  normalized  by  inner  variables.  □,  A,  co3,  Current 
- ,  co3,  KMM;  □,  cO] ,  •,  co3,  Kreplin  and  Eckelmann  (1979). 
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Fig.  15  —  Skewness  of  the  Vorticity  Fluctuations.  D,  w(,  A,  o,  u,,  Current  results. 
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Fig.  16  —  Flatness  of  the  Vorticity  Fluctuations.  □,  u>|  A,  o>2,  o,  Current  results. 
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pig.  17(a)  —  Spanwise  dependence  of  the  correlation  function.  Streamwise  velocity,  a:  x *  =  0.062 
x *  =  125,  Current  results;  — ,  x*  =  5.39,  — ,  x ’  =  149.23,  KMM. 
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Fig.  17(b)  —  Spanwise  dependence  of  the  correlation  function.  Wall-normal  velocity,  a:  2.402,  b:  125 
Current  results;  — ,  jc*  =  5.39, - ,  149.23,  KMM. 
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Fig.  18(b)  —  Streamwise  dependence  of  the  correlation  function.  Wall-normal  velocity,  a:  0.602,  b:  2.402,  c. 
125,  Current  results; . ,  5.39,  — .  149.23,  KMM. 
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Fig.  18(c)  —  Streamwise  dependence  of  the  correlation  function.  Spanwise  velocity,  a:  0.602,  b:  21.06,  c 
125,  Current  results;  — ,  5.39,  — ,  149.23  KMM. 


Fig.  19(a)  -  Wall-normal  dependence  of  the  correlation  coefficient.  Streamwise  velocity,  a:  ,r*  =  0.602.  b 


Fig.  19(b)  —  Wall-normal  dependence  of  the  correlation  coefficient.  Wall-normal  velocity,  a:  0.602,  b:  125. 


normal  dependence  of  the  correlation  function  compared  with  experiments.  Wall-normal 
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Fig.  21  —  Ensemble  average  of  the  streamwise  velocity  for  turbulent  events  detected  with  the  VITA  conditional 
sampling  scheme.  The  detection  parameters  are:  k  =  1.0  and  T*  =  16.  The  detection  is  located  15/*  above 
the  wall. 
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Fig.  25  —  The  normal  velocity  ( u2 )  in  the  vicinity  of  a  low  speed  streak  as  a  function  of  streamwise  location 
and  time,  (x *,/*).  The  figure  is  of  the  velocity  at  a  distance  of  x*  =  15.  The  regions  of  wavelike  disturbance 
are  localized  pockets  of  second  and  fourth  quadrant  Reynolds  stress. 


2700.0 


97 


•  • 


3.0 


Fig.  A2(b)  —  Streamwise  velocity  fluctuation  intensity  profiles.  I:  t  =  0.  2:  t  =  50.  3:  t  =  100.  4:  r 
5:  /  =  200. 
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Fig.  A3(a)  —  Reynolds  number  behavior  in  channel  flow  turbulence  with  and  without  aliasing.  (From  Panda 
(1987))  Wall  Reynolds  number  versus  time  (viscous  units).  33  X  64  x  16  grid:  —  without  dealiasing,  — , 
with  dealiasing. 
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Fig.  A3(b)  —  Reynolds  number  behavior  in  channel  flow  turbulence  with  and  without  aliasing.  (From  Panda 
(1987)).  Centerline  velocity  versus  time. 
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Fig.  B1  —  Root  mean  square  velocity  fluctuation  profiles  computed  from  Equation  Bll. 


